i.atcorr: use same indentation rules for all files

(merge from devbr6, https://trac.osgeo.org/grass/changeset/33216)


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@33217 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.atcorr/6s.cpp

@@ -19,7 +19,7 @@ extern void discom(const GeomCond &geom, const AtmosModel &atms,
 extern void specinterp(const float wl, float& tamoy, float& tamoyp, float& pizmoy, float& pizmoyp,
                        const AerosolConcentration &aerocon, const Altitude &alt);
 extern void enviro (const float difr, const float difa, const float r, const float palt,
-			        const float xmuv, float& fra, float& fae, float& fr);
+		    const float xmuv, float& fra, float& fae, float& fr);
 void printOutput(); // forward declare this function so that it can be used in init_6S
 
 
@@ -38,38 +38,38 @@ static IWave iwave;
 static AtmosModel original_atms;
 int init_6S(char* icnd_name)
 {
-	/* (atmospheric conditions input text file) */
+    /* (atmospheric conditions input text file) */
     ifstream inText;
-	inText.open(icnd_name);
-	if(!inText.is_open()) {
+    inText.open(icnd_name);
+    if(!inText.is_open()) {
         fprintf(stderr, "Unable to open %s\n", icnd_name);
-		return -1;
-	}
+	return -1;
+    }
 
-	/* redirect cin to the input text file */
-	cin.rdbuf(inText.rdbuf());
+    /* redirect cin to the input text file */
+    cin.rdbuf(inText.rdbuf());
 
-	/* read the input geometrical conditions */
-	geom = GeomCond::Parse();
+    /* read the input geometrical conditions */
+    geom = GeomCond::Parse();
 
-	/* read atmospheric model */
-	original_atms = AtmosModel::Parse();
+    /* read atmospheric model */
+    original_atms = AtmosModel::Parse();
     atms = original_atms; /* making a copy */
 
-	/* read aerosol model */
-	aero = AerosolModel::Parse(geom.xmud);
+    /* read aerosol model */
+    aero = AerosolModel::Parse(geom.xmud);
 
     /* read aerosol concentration */
     aerocon = AerosolConcentration::Parse(aero.iaer, atms);
     
-	/* read altitude */
-	alt = Altitude::Parse();
+    /* read altitude */
+    alt = Altitude::Parse();
     alt.init(atms, aerocon);
 
-	/* read iwave stuff */
-	iwave = IWave::Parse();
+    /* read iwave stuff */
+    iwave = IWave::Parse();
    
-	/**********************************************************************c
+    /**********************************************************************c
 	c here, we first compute an equivalent wavelenght which is the input   c
 	c value for monochromatic conditions or the integrated value for a     c
 	c filter function (call equivwl) then, the atmospheric properties are  c
@@ -85,13 +85,13 @@ int init_6S(char* icnd_name)
 	c directional reflectances                                             c
 	c**********************************************************************/
 
-	float wlmoy;
-	if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-	else wlmoy = iwave.wl;
+    float wlmoy;
+    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
+    else wlmoy = iwave.wl;
 
-	discom(geom, atms, aero, aerocon, alt, iwave);
-	float tamoy, tamoyp, pizmoy, pizmoyp;
-	if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
+    discom(geom, atms, aero, aerocon, alt, iwave);
+    float tamoy, tamoyp, pizmoy, pizmoyp;
+    if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
 
     printOutput();
     fflush(stderr);
@@ -107,13 +107,13 @@ void pre_compute_hv(const float height, const float vis)
     alt.set_height(height);
     alt.init(atms, aerocon);
    
-	float wlmoy;
-	if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-	else wlmoy = iwave.wl;
+    float wlmoy;
+    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
+    else wlmoy = iwave.wl;
 
-	discom(geom, atms, aero, aerocon, alt, iwave);
-	float tamoy, tamoyp, pizmoy, pizmoyp;
-	if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
+    discom(geom, atms, aero, aerocon, alt, iwave);
+    float tamoy, tamoyp, pizmoy, pizmoyp;
+    if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
 }
 
 /* Only update those objects that are affected by a visibility change */
@@ -123,13 +123,13 @@ void pre_compute_v(const float vis)
     aerocon.set_visibility(vis, atms);
     alt.init(atms, aerocon);
 
-	float wlmoy;
-	if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-	else wlmoy = iwave.wl;
+    float wlmoy;
+    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
+    else wlmoy = iwave.wl;
 
-	discom(geom, atms, aero, aerocon, alt, iwave);
-	float tamoy, tamoyp, pizmoy, pizmoyp;
-	if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
+    discom(geom, atms, aero, aerocon, alt, iwave);
+    float tamoy, tamoyp, pizmoy, pizmoyp;
+    if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
 }
 
 /* Only update those objects that are affected by a height change */
@@ -139,190 +139,190 @@ void pre_compute_h(const float height)
     alt.set_height(height);
     alt.init(atms, aerocon);
 
-	float wlmoy;
-	if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-	else wlmoy = iwave.wl;
+    float wlmoy;
+    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
+    else wlmoy = iwave.wl;
 
-	discom(geom, atms, aero, aerocon, alt, iwave);
-	float tamoy, tamoyp, pizmoy, pizmoyp;
-	if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
+    discom(geom, atms, aero, aerocon, alt, iwave);
+    float tamoy, tamoyp, pizmoy, pizmoyp;
+    if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
 }
 
 
 void printOutput()
 {
-	static const string head(" 6s version 4.2b ");
+    static const string head(" 6s version 4.2b ");
 
-	cout << endl << endl << endl;
-	Output::Begin(); Output::Repeat(30,'*'); Output::Print(head); Output::Repeat(30,'*'); Output::End();
+    cout << endl << endl << endl;
+    Output::Begin(); Output::Repeat(30,'*'); Output::Print(head); Output::Repeat(30,'*'); Output::End();
 
-	/* ---- geometrical conditions ---- */
-	geom.print();
+    /* ---- geometrical conditions ---- */
+    geom.print();
 
-	/* --- atmospheric model ---- */
-	atms.print();
+    /* --- atmospheric model ---- */
+    atms.print();
 
-	/* --- aerosols model (type) ---- */
-	aero.print();
+    /* --- aerosols model (type) ---- */
+    aero.print();
 
     /* --- aerosols model (concentration) ---- */
     aerocon.print();
 
-	/* --- spectral condition ---- */
-	iwave.print();
+    /* --- spectral condition ---- */
+    iwave.print();
 
-	/* --- ground reflectance (type and spectral variation) ---- */
+    /* --- ground reflectance (type and spectral variation) ---- */
 
-	Output::Ln();
-	Output::WriteLn(22," target type  ");
-	Output::WriteLn(22," -----------  ");
-	Output::WriteLn(10," homogeneous ground ");
+    Output::Ln();
+    Output::WriteLn(22," target type  ");
+    Output::WriteLn(22," -----------  ");
+    Output::WriteLn(10," homogeneous ground ");
 
-	/* 12x a39 f6.3 */
-	static const string reflec[8] = {
-		string(" user defined spectral reflectance     "),
-		string(" monochromatic reflectance "),
-		string(" constant reflectance over the spectra "),
-		string(" spectral vegetation ground reflectance"),
-		string(" spectral clear water reflectance      "),
-		string(" spectral dry sand ground reflectance  "),
-		string(" spectral lake water reflectance       "),
-		string(" spectral volcanic debris reflectance  ")
-	};
+    /* 12x a39 f6.3 */
+    static const string reflec[8] = {
+	string(" user defined spectral reflectance     "),
+	string(" monochromatic reflectance "),
+	string(" constant reflectance over the spectra "),
+	string(" spectral vegetation ground reflectance"),
+	string(" spectral clear water reflectance      "),
+	string(" spectral dry sand ground reflectance  "),
+	string(" spectral lake water reflectance       "),
+	string(" spectral volcanic debris reflectance  ")
+    };
 
     float rocave = 0;       /* block of code in Fortran will always compute 0 */
+    ostringstream s;
+    s.setf(ios::fixed, ios::floatfield);
+    s << setprecision(3);
+    s << reflec[2] << setw(9) << rocave << ends;
+    Output::WriteLn(12, s.str());
+
+
+    /* --- pressure at ground level (174) and altitude (175) ---- */
+    Output::Ln();
+    Output::WriteLn(22," target elevation description ");
+    Output::WriteLn(22," ---------------------------- ");
+
+    ostringstream s1;
+    s1.setf(ios::fixed, ios::floatfield);
+    s1 << setprecision(2);
+    s1 << " ground pressure  [mb]     " << setw(9) << atms.p[0] << ends;
+    Output::WriteLn(10,s1.str());
+
+    ostringstream s2;
+    s2.setf(ios::fixed, ios::floatfield);
+    s2 << setprecision(3);
+    s2 << " ground altitude  [km]    " << setw(9) << alt.xps << ends;
+    Output::WriteLn(10,s2.str());
+
+    if( alt.xps > 0 )
+    {
+	Output::WriteLn(15," gaseous content at target level: ");
+
 	ostringstream s;
 	s.setf(ios::fixed, ios::floatfield);
 	s << setprecision(3);
-	s << reflec[2] << setw(9) << rocave << ends;
-	Output::WriteLn(12, s.str());
-
-
-	/* --- pressure at ground level (174) and altitude (175) ---- */
-	Output::Ln();
-	Output::WriteLn(22," target elevation description ");
-	Output::WriteLn(22," ---------------------------- ");
-
-	ostringstream s1;
-	s1.setf(ios::fixed, ios::floatfield);
-	s1 << setprecision(2);
-	s1 << " ground pressure  [mb]     " << setw(9) << atms.p[0] << ends;
-	Output::WriteLn(10,s1.str());
-
-	ostringstream s2;
-	s2.setf(ios::fixed, ios::floatfield);
-	s2 << setprecision(3);
-	s2 << " ground altitude  [km]    " << setw(9) << alt.xps << ends;
-	Output::WriteLn(10,s2.str());
-
-	if( alt.xps > 0 )
-	{
-		Output::WriteLn(15," gaseous content at target level: ");
-
-		ostringstream s;
-		s.setf(ios::fixed, ios::floatfield);
-		s << setprecision(3);
-		s << " uh2o=" << setw(9) << atms.uw << " g/cm2      "
-		  << "  uo3=" << setw(9) << atms.uo3 << " cm-atm" << ends;
-		Output::WriteLn(15,s.str());
-	}
+	s << " uh2o=" << setw(9) << atms.uw << " g/cm2      "
+	  << "  uo3=" << setw(9) << atms.uo3 << " cm-atm" << ends;
+	Output::WriteLn(15,s.str());
+    }
 
-	alt.print();
+    alt.print();
 
-	/* ---- atmospheric correction  ---- */
-	Output::Ln();
-	Output::WriteLn(23," atmospheric correction activated ");
-	Output::WriteLn(23," -------------------------------- ");
+    /* ---- atmospheric correction  ---- */
+    Output::Ln();
+    Output::WriteLn(23," atmospheric correction activated ");
+    Output::WriteLn(23," -------------------------------- ");
 }
 
 TransformInput compute()
 {
-	const float accu3 = 1e-07;
+    const float accu3 = 1e-07;
 /* ---- initilialization	 very liberal :) */
-	int i, j;
-
-	float fr = 0;
-	float rad = 0;
-	float sb = 0;
-	float seb = 0;
-	float refet = 0;
-	float refet1 = 0;
-	float refet2 = 0;
-	float refet3 = 0;
-	float alumet = 0;
-	float tgasm = 0;
-	float rog = 0;
-	float dgasm = 0;
-	float ugasm = 0;
-	float sdwava = 0;
-	float sdozon = 0;
-	float sddica = 0;
-	float sdoxyg = 0;
-	float sdniox = 0;
-	float sdmoca = 0;
-	float sdmeth = 0;
-
-	float suwava = 0;
-	float suozon = 0;
-	float sudica = 0;
-	float suoxyg = 0;
-	float suniox = 0;
-	float sumoca = 0;
-	float sumeth = 0;
-	float stwava = 0;
-	float stozon = 0;
-	float stdica = 0;
-	float stoxyg = 0;
-	float stniox = 0;
-	float stmoca = 0;
-	float stmeth = 0;
-	float sodray = 0;
-	float sodrayp = 0;
-	float sodaer = 0;
-	float sodaerp = 0;
-	float sodtot = 0;
-	float sodtotp = 0;
-	float fophsr = 0;
-	float fophsa = 0;
-	float sroray = 0;
-	float sroaer = 0;
-	float srotot = 0;
-	float ssdaer = 0;
-	float sdtotr = 0;
-	float sdtota = 0;
-	float sdtott = 0;
-	float sutotr = 0;
-	float sutota = 0;
-	float sutott = 0;
-	float sasr = 0;
-	float sasa = 0;
-	float sast = 0;
-
-	float ani[2][3];
-	float aini[2][3];
-	float anr[2][3];
-	float ainr[2][3];
-
-	for(i = 0; i < 2; i++)
-		for(j = 0; j < 3; j++)
-		{
-			ani[i][j] = 0;
-			aini[i][j] = 0;
-			anr[i][j] = 0;
-			ainr[i][j] = 0;
-		}
-
-	/* ---- spectral loop ---- */
-    if (iwave.iwave == -2)
+    int i, j;
+
+    float fr = 0;
+    float rad = 0;
+    float sb = 0;
+    float seb = 0;
+    float refet = 0;
+    float refet1 = 0;
+    float refet2 = 0;
+    float refet3 = 0;
+    float alumet = 0;
+    float tgasm = 0;
+    float rog = 0;
+    float dgasm = 0;
+    float ugasm = 0;
+    float sdwava = 0;
+    float sdozon = 0;
+    float sddica = 0;
+    float sdoxyg = 0;
+    float sdniox = 0;
+    float sdmoca = 0;
+    float sdmeth = 0;
+
+    float suwava = 0;
+    float suozon = 0;
+    float sudica = 0;
+    float suoxyg = 0;
+    float suniox = 0;
+    float sumoca = 0;
+    float sumeth = 0;
+    float stwava = 0;
+    float stozon = 0;
+    float stdica = 0;
+    float stoxyg = 0;
+    float stniox = 0;
+    float stmoca = 0;
+    float stmeth = 0;
+    float sodray = 0;
+    float sodrayp = 0;
+    float sodaer = 0;
+    float sodaerp = 0;
+    float sodtot = 0;
+    float sodtotp = 0;
+    float fophsr = 0;
+    float fophsa = 0;
+    float sroray = 0;
+    float sroaer = 0;
+    float srotot = 0;
+    float ssdaer = 0;
+    float sdtotr = 0;
+    float sdtota = 0;
+    float sdtott = 0;
+    float sutotr = 0;
+    float sutota = 0;
+    float sutott = 0;
+    float sasr = 0;
+    float sasa = 0;
+    float sast = 0;
+
+    float ani[2][3];
+    float aini[2][3];
+    float anr[2][3];
+    float ainr[2][3];
+
+    for(i = 0; i < 2; i++)
+	for(j = 0; j < 3; j++)
 	{
-		Output::WriteLn(1,"wave   total  total  total  total  atm.   swl    step   sbor   dsol   toar ");
-		Output::WriteLn(1,"       gas    scat   scat   spheri intr   ");
-		Output::WriteLn(1,"       trans  down   up     albedo refl   ");
+	    ani[i][j] = 0;
+	    aini[i][j] = 0;
+	    anr[i][j] = 0;
+	    ainr[i][j] = 0;
 	}
 
-	int l;
-	for(l = iwave.iinf; l <= iwave.isup; l++)
-	{
+    /* ---- spectral loop ---- */
+    if (iwave.iwave == -2)
+    {
+	Output::WriteLn(1,"wave   total  total  total  total  atm.   swl    step   sbor   dsol   toar ");
+	Output::WriteLn(1,"       gas    scat   scat   spheri intr   ");
+	Output::WriteLn(1,"       trans  down   up     albedo refl   ");
+    }
+
+    int l;
+    for(l = iwave.iinf; l <= iwave.isup; l++)
+    {
         float sbor = iwave.ffu.s[l];
 
         if(l == iwave.iinf || l == iwave.isup) sbor *= 0.5f;
@@ -332,16 +332,16 @@ TransformInput compute()
         float roe = 0; /* roel[l]; */
         float wl = 0.25f + l * step;
 
-		AbstraStruct as;
-		float uwus, uo3us;		/* initialized in abstra */
+	AbstraStruct as;
+	float uwus, uo3us;		/* initialized in abstra */
 
-		abstra(atms, alt, wl, (float)geom.xmus, (float)geom.xmuv, atms.uw / 2.0f, atms.uo3,
-			   uwus, uo3us, alt.puw / 2.0f, alt.puo3, alt.puwus, alt.puo3us, as);
+	abstra(atms, alt, wl, (float)geom.xmus, (float)geom.xmuv, atms.uw / 2.0f, atms.uo3,
+	       uwus, uo3us, alt.puw / 2.0f, alt.puo3, alt.puwus, alt.puo3us, as);
 
-		float attwava = as.ttwava;
+	float attwava = as.ttwava;
 
-		abstra(atms, alt, wl, (float)geom.xmus, (float)geom.xmuv, atms.uw, atms.uo3,
-			   uwus, uo3us, alt.puw, alt.puo3, alt.puwus, alt.puo3us, as);
+	abstra(atms, alt, wl, (float)geom.xmus, (float)geom.xmuv, atms.uw, atms.uo3,
+	       uwus, uo3us, alt.puw, alt.puo3, alt.puwus, alt.puo3us, as);
 
         if (as.dtwava < accu3) as.dtwava = 0;
         if (as.dtozon < accu3) as.dtozon = 0;
@@ -366,9 +366,9 @@ TransformInput compute()
         swl = swl * geom.dsol;
         float coef = sbor * step * swl;
 
-		InterpStruct is;
-		memset(&is, 0, sizeof(is));
-		interp(aero.iaer, alt.idatmp, wl, aerocon.taer55, alt.taer55p, (float)geom.xmud, is);
+	InterpStruct is;
+	memset(&is, 0, sizeof(is));
+	interp(aero.iaer, alt.idatmp, wl, aerocon.taer55, alt.taer55p, (float)geom.xmud, is);
 
 
         float dgtot = as.dtwava * as.dtozon * as.dtdica * as.dtoxyg * as.dtniox * as.dtmeth * as.dtmoca;
@@ -380,12 +380,12 @@ TransformInput compute()
         float edifa = (float)(is.utota - exp(-is.taerp / geom.xmuv));
 
 
-		float fra, fae;
-		enviro(edifr, edifa, rad, alt.palt, (float)geom.xmuv, fra, fae, fr);
+	float fra, fae;
+	enviro(edifr, edifa, rad, alt.palt, (float)geom.xmuv, fra, fae, fr);
 
-		float avr = roc * fr + (1 - fr) * roe;
-		float rsurf = (float)(roc * is.dtott * exp(-(is.trayp + is.taerp) / geom.xmuv) / (1 - avr * is.astot)
-					+ avr * is.dtott * (is.utott - exp(-(is.trayp + is.taerp) / geom.xmuv)) / (1 - avr * is.astot));
+	float avr = roc * fr + (1 - fr) * roe;
+	float rsurf = (float)(roc * is.dtott * exp(-(is.trayp + is.taerp) / geom.xmuv) / (1 - avr * is.astot)
+			      + avr * is.dtott * (is.utott - exp(-(is.trayp + is.taerp) / geom.xmuv)) / (1 - avr * is.astot));
         float ratm1 = (is.romix - is.rorayl) * tgtot + is.rorayl * tgp1;
         float ratm3 = is.romix * tgp1;
         float ratm2 = (is.romix - is.rorayl) * tgp2 + is.rorayl * tgp1;
@@ -393,24 +393,24 @@ TransformInput compute()
         float romeas2 = ratm2 + rsurf * tgtot;
         float romeas3 = ratm3 + rsurf * tgtot;
 
-		/* computing integrated values over the spectral band */
+	/* computing integrated values over the spectral band */
         if (iwave.iwave == -2)
-		{
-			Output::Begin();
-			ostringstream s;
-			s.setf(ios::fixed, ios::floatfield);
-			s.precision(4);
-			s	<< setw(10) << wl << " "
-				<< setw(10) << tgtot << " "
-				<< setw(10) << is.dtott << " "
-				<< setw(10) << is.utott << " "
-				<< setw(10) << is.astot << " "
-				<< setw(10) << ratm2 << " "
-				<< setprecision(1) << setw(7) << swl << " "
-				<< setprecision(4) << setw(10) << step << " "
-				<< setw(10) << sbor << " "
-				<< setw(10) << geom.dsol << " "
-				<< setw(10) << romeas2;
+	{
+	    Output::Begin();
+	    ostringstream s;
+	    s.setf(ios::fixed, ios::floatfield);
+	    s.precision(4);
+	    s	<< setw(10) << wl << " "
+		<< setw(10) << tgtot << " "
+		<< setw(10) << is.dtott << " "
+		<< setw(10) << is.utott << " "
+		<< setw(10) << is.astot << " "
+		<< setw(10) << ratm2 << " "
+		<< setprecision(1) << setw(7) << swl << " "
+		<< setprecision(4) << setw(10) << step << " "
+		<< setw(10) << sbor << " "
+		<< setw(10) << geom.dsol << " "
+		<< setw(10) << romeas2;
         }
 
         
@@ -469,7 +469,7 @@ TransformInput compute()
         sb = sb + sbor * step;
         seb = seb + coef;
 
-		/* output at the ground level. */
+	/* output at the ground level. */
         float tdir = (float)exp(-(is.tray + is.taer) / geom.xmus);
         float tdif = is.dtott - tdir;
         float etn = is.dtott * dgtot / (1 - avr * is.astot);
@@ -481,16 +481,16 @@ TransformInput compute()
         float ee0 = (float)(geom.xmus * swl * dgtot * avr * is.astot * is.dtott / (1 - avr * is.astot));
 
         if (etn > accu3)
-		{
-           ani[0][0] = esn / etn;
-           ani[0][1] = ea0n / etn;
-           ani[0][2] = ee0n / etn;
-		}
+	{
+	    ani[0][0] = esn / etn;
+	    ani[0][1] = ea0n / etn;
+	    ani[0][2] = ee0n / etn;
+	}
         else
-		{
-           ani[0][0] = 0;
-           ani[0][1] = 0;
-           ani[0][2] = 0;
+	{
+	    ani[0][0] = 0;
+	    ani[0][1] = 0;
+	    ani[0][2] = 0;
         }
 
         ani[1][0] = es;
@@ -498,13 +498,13 @@ TransformInput compute()
         ani[1][2] = ee0;
 
 
-		for(j = 0; j < 3; j++)
-		{
-          aini[0][j] = aini[0][j] + ani[0][j] * coef;
-          aini[1][j] = aini[1][j] + ani[1][j] * sbor * step;
-		}
+	for(j = 0; j < 3; j++)
+	{
+	    aini[0][j] = aini[0][j] + ani[0][j] * coef;
+	    aini[1][j] = aini[1][j] + ani[1][j] * sbor * step;
+	}
 
-		/* output at satellite level */
+	/* output at satellite level */
         float tmdir = (float)exp(-(is.tray + is.taer) / geom.xmuv);
         float tmdif = is.utott - tmdir;
         float xla0n = ratm2;
@@ -520,78 +520,78 @@ TransformInput compute()
         anr[1][1] = xle;
         anr[1][2] = xlt;
 
-		for(j = 0; j < 3; j++)
-		{
-          ainr[0][j] = ainr[0][j] + anr[0][j] * coef;
-          ainr[1][j] = ainr[1][j] + anr[1][j] * sbor * step;
-		}
+	for(j = 0; j < 3; j++)
+	{
+	    ainr[0][j] = ainr[0][j] + anr[0][j] * coef;
+	    ainr[1][j] = ainr[1][j] + anr[1][j] * sbor * step;
 	}
+    }
 
 
-	/* ---- integrated values of apparent reflectance, radiance          ----*/
-	/* ---- and gaseous transmittances (total,downward,separately gases) ----*/
+    /* ---- integrated values of apparent reflectance, radiance          ----*/
+    /* ---- and gaseous transmittances (total,downward,separately gases) ----*/
     refet = refet / seb;
-	refet1 = refet1 / seb;
-	refet2 = refet2 / seb;
-	refet3 = refet3 / seb;
-	tgasm = tgasm / seb;
-	dgasm = dgasm / seb;
-	ugasm = ugasm / seb;
-	sasa = sasa / seb;
-	sasr = sasr / seb;
-	sast = sast / seb;
-	sdniox = sdniox / seb;
-	sdmoca = sdmoca / seb;
-	sdmeth = sdmeth / seb;
-	sdwava = sdwava / seb;
-	sdozon = sdozon / seb;
-	sddica = sddica / seb;
-	suniox = suniox / seb;
-	sumoca = sumoca / seb;
-	sumeth = sumeth / seb;
-	suwava = suwava / seb;
-	suozon = suozon / seb;
-	sudica = sudica / seb;
-	suoxyg = suoxyg / seb;
-	sdoxyg = sdoxyg / seb;
-	stniox = stniox / seb;
-	stmoca = stmoca / seb;
-	stmeth = stmeth / seb;
-	stwava = stwava / seb;
-	stozon = stozon / seb;
-	stdica = stdica / seb;
-	stoxyg = stoxyg / seb;
-	sdtotr = sdtotr / seb;
-	sdtota = sdtota / seb;
-	sdtott = sdtott / seb;
-
-	sutotr = sutotr / seb;
-	sutota = sutota / seb;
-	sutott = sutott / seb;
-	rog = rog / seb;
-	sroray = sroray / seb;
-	sroaer = sroaer / seb;
-	srotot = srotot / seb;
-	alumet = alumet / sb;
-	float pizera = 0.0f;
-	if(aero.iaer != 0) pizera = ssdaer / sodaer;
-	sodray = sodray / seb;
-	sodaer = sodaer / seb;
-	sodtot = sodtot / seb;
-	sodrayp = sodrayp / seb;
-	sodaerp = sodaerp / seb;
-	sodtotp = sodtotp / seb;
-	fophsa = fophsa / seb;
-	fophsr = fophsr / seb;
-
-	for(j = 0; j < 3; j++)
-
-	{
+    refet1 = refet1 / seb;
+    refet2 = refet2 / seb;
+    refet3 = refet3 / seb;
+    tgasm = tgasm / seb;
+    dgasm = dgasm / seb;
+    ugasm = ugasm / seb;
+    sasa = sasa / seb;
+    sasr = sasr / seb;
+    sast = sast / seb;
+    sdniox = sdniox / seb;
+    sdmoca = sdmoca / seb;
+    sdmeth = sdmeth / seb;
+    sdwava = sdwava / seb;
+    sdozon = sdozon / seb;
+    sddica = sddica / seb;
+    suniox = suniox / seb;
+    sumoca = sumoca / seb;
+    sumeth = sumeth / seb;
+    suwava = suwava / seb;
+    suozon = suozon / seb;
+    sudica = sudica / seb;
+    suoxyg = suoxyg / seb;
+    sdoxyg = sdoxyg / seb;
+    stniox = stniox / seb;
+    stmoca = stmoca / seb;
+    stmeth = stmeth / seb;
+    stwava = stwava / seb;
+    stozon = stozon / seb;
+    stdica = stdica / seb;
+    stoxyg = stoxyg / seb;
+    sdtotr = sdtotr / seb;
+    sdtota = sdtota / seb;
+    sdtott = sdtott / seb;
+
+    sutotr = sutotr / seb;
+    sutota = sutota / seb;
+    sutott = sutott / seb;
+    rog = rog / seb;
+    sroray = sroray / seb;
+    sroaer = sroaer / seb;
+    srotot = srotot / seb;
+    alumet = alumet / sb;
+    float pizera = 0.0f;
+    if(aero.iaer != 0) pizera = ssdaer / sodaer;
+    sodray = sodray / seb;
+    sodaer = sodaer / seb;
+    sodtot = sodtot / seb;
+    sodrayp = sodrayp / seb;
+    sodaerp = sodaerp / seb;
+    sodtotp = sodtotp / seb;
+    fophsa = fophsa / seb;
+    fophsr = fophsr / seb;
+
+    for(j = 0; j < 3; j++)
+
+    {
         aini[0][j] = aini[0][j] / seb;
         ainr[0][j] = ainr[0][j] / seb;
         aini[1][j] = aini[1][j] / sb;
         ainr[1][j] = ainr[1][j] / sb;
-	}
+    }
 
     /* Prepare data for final dn transformation */
     TransformInput ti;

imagery/i.atcorr/AerosolConcentration.cpp

@@ -24,119 +24,119 @@ void AerosolConcentration::parse(const long int _iaer, const AtmosModel& atms)
 {
     iaer = _iaer;
 
-	taer55 = 0.f;
-	cin >> v;
-	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
+    taer55 = 0.f;
+    cin >> v;
+    cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
 
-	if(v == 0)
-	{
-		cin >> taer55;
-		cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
-		v = (float)(exp(-log(taer55/2.7628f)/0.79902f));
-	}
-	else if(v > 0) oda550(v, atms);
+    if(v == 0)
+    {
+	cin >> taer55;
+	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
+	v = (float)(exp(-log(taer55/2.7628f)/0.79902f));
+    }
+    else if(v > 0) oda550(v, atms);
 }
 
 void AerosolConcentration::oda550(const float v, const AtmosModel& atms)
 {
-	/* aerosol optical depth at wl=550nm */
-	/* vertical repartition of aerosol density for v=23km */
-	/*                ( in nbr of part/cm3 ) */
-	static const float an23[34] = {
-		2.828e+03,1.244e+03,5.371e+02,2.256e+02,1.192e+02,
-		8.987e+01,6.337e+01,5.890e+01,6.069e+01,5.818e+01,
-		5.675e+01,5.317e+01,5.585e+01,5.156e+01,5.048e+01,
-		4.744e+01,4.511e+01,4.458e+01,4.314e+01,3.634e+01,
-		2.667e+01,1.933e+01,1.455e+01,1.113e+01,8.826e+00,
-		7.429e+00,2.238e+00,5.890e-01,1.550e-01,4.082e-02,
-		1.078e-02,5.550e-05,1.969e-08,0.000e+00
-	};
-
-
-	/* vertical repartition of aerosol density for v=5km */
-	/*                ( in nbr of part/cm3 ) */
-	static const float an5[34] = {
-		1.378e+04,5.030e+03,1.844e+03,6.731e+02,2.453e+02,
-		8.987e+01,6.337e+01,5.890e+01,6.069e+01,5.818e+01,
-		5.675e+01,5.317e+01,5.585e+01,5.156e+01,5.048e+01,
-		4.744e+01,4.511e+01,4.458e+01,4.314e+01,3.634e+01,
-		2.667e+01,1.933e+01,1.455e+01,1.113e+01,8.826e+00,
-		7.429e+00,2.238e+00,5.890e-01,1.550e-01,4.082e-02,
-		1.078e-02,5.550e-05,1.969e-08,0.000e+00
-	};
-
-	taer55 = 0;
-	if(fabs(v) <= 0) return;
-	if(iaer == 0) return;
-
-	for(int k = 0; k < 32; k++)
-	{
-		float dz = atms.z[k+1] - atms.z[k];
-		float az = (115.f / 18.f) * (an5[k] - an23[k]);
-		float az1 = (115.f / 18.f) * (an5[k+1] - an23[k+1]);
-
-		float bz = (5.f * an5[k] / 18.f) - (23.f * an23[k] / 18.f);
-		float bz1 = (5.f * an5[k+1] / 18.f) - (23.f * an23[k+1] / 18.f);
-
-		float bnz = az / v - bz;
-		float bnz1 = az1 / v - bz1;
-
-		float ev = (float)(dz * exp((log(bnz) + log(bnz1)) / 2));
-		taer55 += ev * sigma * 1.0e-03f;
-	}
+    /* aerosol optical depth at wl=550nm */
+    /* vertical repartition of aerosol density for v=23km */
+    /*                ( in nbr of part/cm3 ) */
+    static const float an23[34] = {
+	2.828e+03,1.244e+03,5.371e+02,2.256e+02,1.192e+02,
+	8.987e+01,6.337e+01,5.890e+01,6.069e+01,5.818e+01,
+	5.675e+01,5.317e+01,5.585e+01,5.156e+01,5.048e+01,
+	4.744e+01,4.511e+01,4.458e+01,4.314e+01,3.634e+01,
+	2.667e+01,1.933e+01,1.455e+01,1.113e+01,8.826e+00,
+	7.429e+00,2.238e+00,5.890e-01,1.550e-01,4.082e-02,
+	1.078e-02,5.550e-05,1.969e-08,0.000e+00
+    };
+
+
+    /* vertical repartition of aerosol density for v=5km */
+    /*                ( in nbr of part/cm3 ) */
+    static const float an5[34] = {
+	1.378e+04,5.030e+03,1.844e+03,6.731e+02,2.453e+02,
+	8.987e+01,6.337e+01,5.890e+01,6.069e+01,5.818e+01,
+	5.675e+01,5.317e+01,5.585e+01,5.156e+01,5.048e+01,
+	4.744e+01,4.511e+01,4.458e+01,4.314e+01,3.634e+01,
+	2.667e+01,1.933e+01,1.455e+01,1.113e+01,8.826e+00,
+	7.429e+00,2.238e+00,5.890e-01,1.550e-01,4.082e-02,
+	1.078e-02,5.550e-05,1.969e-08,0.000e+00
+    };
+
+    taer55 = 0;
+    if(fabs(v) <= 0) return;
+    if(iaer == 0) return;
+
+    for(int k = 0; k < 32; k++)
+    {
+	float dz = atms.z[k+1] - atms.z[k];
+	float az = (115.f / 18.f) * (an5[k] - an23[k]);
+	float az1 = (115.f / 18.f) * (an5[k+1] - an23[k+1]);
+
+	float bz = (5.f * an5[k] / 18.f) - (23.f * an23[k] / 18.f);
+	float bz1 = (5.f * an5[k+1] / 18.f) - (23.f * an23[k+1] / 18.f);
+
+	float bnz = az / v - bz;
+	float bnz1 = az1 / v - bz1;
+
+	float ev = (float)(dz * exp((log(bnz) + log(bnz1)) / 2));
+	taer55 += ev * sigma * 1.0e-03f;
+    }
 }
 
 void AerosolConcentration::print()
 {
-	/* --- aerosol model (concentration) ---- */
+    /* --- aerosol model (concentration) ---- */
+    Output::Begin();
+    Output::End();
+    if(iaer == 0) return;
+
+    Output::Begin();
+    Output::Repeat(10, ' ');
+    Output::Print(" optical condition identity :");
+    Output::End();
+    if(fabs(v) <= xacc)
+    {
 	Output::Begin();
+	Output::Repeat(15, ' ');
+	Output::Print(" user def. opt. thick. at 550nm :");
+	ostringstream s;
+	s.setf(ios::fixed, ios::floatfield);
+	s << setprecision(4);
+	s << setw(11) << taer55 << ends;
+	Output::Print(s.str());
 	Output::End();
-	if(iaer == 0) return;
-
+    }
+    else
+    {
 	Output::Begin();
-	Output::Repeat(10, ' ');
-	Output::Print(" optical condition identity :");
+	Output::Repeat(15, ' ');
+	Output::Print(" visibility :");
+	ostringstream s1;
+	s1.setf(ios::fixed, ios::floatfield);
+	s1 << setprecision(2);
+	s1 << setw(8) << v << ends;
+	Output::Print(s1.str());
+	Output::Print(" km  opt. thick. 550nm :");
+	ostringstream s2;
+	s2.setf(ios::fixed, ios::floatfield);
+	s2 << setprecision(4);
+	s2 << setw(9) << taer55 << ends;
+	Output::Print(s2.str());
 	Output::End();
-	if(fabs(v) <= xacc)
-	{
-		Output::Begin();
-		Output::Repeat(15, ' ');
-		Output::Print(" user def. opt. thick. at 550nm :");
-		ostringstream s;
-		s.setf(ios::fixed, ios::floatfield);
-		s << setprecision(4);
-		s << setw(11) << taer55 << ends;
-		Output::Print(s.str());
-		Output::End();
-	}
-	else
-	{
-		Output::Begin();
-		Output::Repeat(15, ' ');
-		Output::Print(" visibility :");
-		ostringstream s1;
-		s1.setf(ios::fixed, ios::floatfield);
-		s1 << setprecision(2);
-		s1 << setw(8) << v << ends;
-		Output::Print(s1.str());
-		Output::Print(" km  opt. thick. 550nm :");
-		ostringstream s2;
-		s2.setf(ios::fixed, ios::floatfield);
-		s2 << setprecision(4);
-		s2 << setw(9) << taer55 << ends;
-		Output::Print(s2.str());
-		Output::End();
-	}
+    }
 
-	Output::Begin();
-	Output::End();    
+    Output::Begin();
+    Output::End();    
 }
 
 
 AerosolConcentration AerosolConcentration::Parse(const long int iaer, const AtmosModel& atms)
 {
-	AerosolConcentration aerocon;
-	aerocon.parse(iaer, atms);
-	return aerocon;
+    AerosolConcentration aerocon;
+    aerocon.parse(iaer, atms);
+    return aerocon;
 }
 

imagery/i.atcorr/Altitude.cpp

@@ -7,186 +7,186 @@
    is not at sea level.
 
 
-	Given the altitude of the target in kilometers as input, we transform the
-	original atmospheric profile (Pressure, Temperature, Water Vapor, Ozone) 
-	so that first level of the new profile is the one at the target altitude. 
-	We also compute the new integrated content in water vapor and ozone, that
-	are used as outputs or in computations when the user chooses to enter a
-	specific amount of Ozone and Water Vapor.
+   Given the altitude of the target in kilometers as input, we transform the
+   original atmospheric profile (Pressure, Temperature, Water Vapor, Ozone) 
+   so that first level of the new profile is the one at the target altitude. 
+   We also compute the new integrated content in water vapor and ozone, that
+   are used as outputs or in computations when the user chooses to enter a
+   specific amount of Ozone and Water Vapor.
 */
 void Altitude::pressure(AtmosModel& atms, float& uw, float& uo3)
 {
-	/* log linear interpolation */
-	if(xps >= 100) xps = 99.99f;
+    /* log linear interpolation */
+    if(xps >= 100) xps = 99.99f;
 		
-	int i;
-	for(i = 0; atms.z[i] <= xps; i++);
+    int i;
+    for(i = 0; atms.z[i] <= xps; i++);
 		
-	int isup = i;
-	int iinf = i - 1;
-
-	float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
-	float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
-	float ps = (float)exp((xps - xb) / xa);
-
-	/* interpolating temperature wator vapor and ozone profile versus altitude */
-	float xalt = xps;
-	float xtemp = (atms.t[isup] - atms.t[iinf]) / (atms.z[isup] - atms.z[iinf]);
-	xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
-	float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
-	xwo = xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
-	float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
-	xwh = xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];
-
-	/* updating atmospheric profile
-	  1rst level: target     , complete to 34
-	  with interpolated layers */
-	atms.z[0] = xalt;                                                          
-	atms.p[0] = ps;
-	atms.t[0] = xtemp;
-	atms.wh[0] = xwh;
-	atms.wo[0] = xwo;
-
-	for (i = 1; i < 33 - iinf; ++i)
+    int isup = i;
+    int iinf = i - 1;
+
+    float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
+    float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
+    float ps = (float)exp((xps - xb) / xa);
+
+    /* interpolating temperature wator vapor and ozone profile versus altitude */
+    float xalt = xps;
+    float xtemp = (atms.t[isup] - atms.t[iinf]) / (atms.z[isup] - atms.z[iinf]);
+    xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
+    float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
+    xwo = xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
+    float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
+    xwh = xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];
+
+    /* updating atmospheric profile
+       1rst level: target     , complete to 34
+       with interpolated layers */
+    atms.z[0] = xalt;                                                          
+    atms.p[0] = ps;
+    atms.t[0] = xtemp;
+    atms.wh[0] = xwh;
+    atms.wo[0] = xwo;
+
+    for (i = 1; i < 33 - iinf; ++i)
     {
-		atms.z[i] = atms.z[i + iinf];
-		atms.p[i] = atms.p[i + iinf];
-		atms.t[i] = atms.t[i + iinf];
-		atms.wh[i] = atms.wh[i + iinf];
-		atms.wo[i] = atms.wo[i + iinf];
+	atms.z[i] = atms.z[i + iinf];
+	atms.p[i] = atms.p[i + iinf];
+	atms.t[i] = atms.t[i + iinf];
+	atms.wh[i] = atms.wh[i + iinf];
+	atms.wo[i] = atms.wo[i + iinf];
     }
-	int l = 33 - iinf - 1;
-	for (i = l; i < 34; ++i)
+    int l = 33 - iinf - 1;
+    for (i = l; i < 34; ++i)
     {
-      atms.z[i] = (atms.z[33] - atms.z[l]) * (i - l) / (33 - l) + atms.z[l];
-      atms.p[i] = (atms.p[33] - atms.p[l]) * (i - l) / (33 - l) + atms.p[l];
-      atms.t[i] = (atms.t[33] - atms.t[l]) * (i - l) / (33 - l) + atms.t[l];
-      atms.wh[i] = (atms.wh[33] - atms.wh[l]) * (i - l) / (33 - l) + atms.wh[l];
-      atms.wo[i] = (atms.wo[33] - atms.wo[l]) * (i - l) / (33 - l) + atms.wo[l];
+	atms.z[i] = (atms.z[33] - atms.z[l]) * (i - l) / (33 - l) + atms.z[l];
+	atms.p[i] = (atms.p[33] - atms.p[l]) * (i - l) / (33 - l) + atms.p[l];
+	atms.t[i] = (atms.t[33] - atms.t[l]) * (i - l) / (33 - l) + atms.t[l];
+	atms.wh[i] = (atms.wh[33] - atms.wh[l]) * (i - l) / (33 - l) + atms.wh[l];
+	atms.wo[i] = (atms.wo[33] - atms.wo[l]) * (i - l) / (33 - l) + atms.wo[l];
     }
 
-	/* compute modified h2o and o3 integrated content */
-	uw = 0;
-	uo3 = 0;
-	const float g = 98.1f;
-	const float air = 0.028964f/0.0224f;
-	const float ro3 = 0.048f/0.0224f;
-
-	float rmwh[34];
-	float rmo3[34];
-	int k;
-	for (k = 0; k < 33; ++k)
+    /* compute modified h2o and o3 integrated content */
+    uw = 0;
+    uo3 = 0;
+    const float g = 98.1f;
+    const float air = 0.028964f/0.0224f;
+    const float ro3 = 0.048f/0.0224f;
+
+    float rmwh[34];
+    float rmo3[34];
+    int k;
+    for (k = 0; k < 33; ++k)
     {
-		float roair = air * 273.16f * atms.p[k] / (atms.t[k] * 1013.25f);
-		rmwh[k] = atms.wh[k] / (roair * 1e3f);
-		rmo3[k] = atms.wo[k] / (roair * 1e3f);
+	float roair = air * 273.16f * atms.p[k] / (atms.t[k] * 1013.25f);
+	rmwh[k] = atms.wh[k] / (roair * 1e3f);
+	rmo3[k] = atms.wo[k] / (roair * 1e3f);
     }
 
-	for (k = 1; k < 33; ++k)
+    for (k = 1; k < 33; ++k)
     {
-		float ds = (atms.p[k - 1] - atms.p[k]) / atms.p[0];
-		uw += (rmwh[k] + rmwh[k - 1]) * ds / 2.f;
-		uo3 += (rmo3[k] + rmo3[k - 1]) * ds / 2.f;
+	float ds = (atms.p[k - 1] - atms.p[k]) / atms.p[0];
+	uw += (rmwh[k] + rmwh[k - 1]) * ds / 2.f;
+	uo3 += (rmo3[k] + rmo3[k - 1]) * ds / 2.f;
     }
-	uw = uw * atms.p[0] * 100.f / g;
-	uo3 = uo3 * atms.p[0] * 100.f / g;
-	uo3 = uo3 * 1e3f / ro3;
+    uw = uw * atms.p[0] * 100.f / g;
+    uo3 = uo3 * atms.p[0] * 100.f / g;
+    uo3 = uo3 * 1e3f / ro3;
 }
 
 /*
-Function: Update the atmospheric profile (P(z),T(z),H2O(z),O3(z)) in case the observer is on
-board an aircraft.
-
-Description: Given the altitude or pressure at aircraft level as input, the first task is to
-compute the altitude (in case the pressure has been entered) or the pressure (in case the altitude has
-been entered) at plane level. Then, a new atmospheric profile is created (Pp,Tp,H2Op,O3p) for which
-the last level is located at the plane altitude. This profile is used in the gaseous absorption
-computation (ABSTRA.f) for the path from target to sensor (upward transmission). The ozone and
-water vapor integrated content of the "plane" atmospheric profile are also an output of this
-subroutine. The last output is the proportion of molecules below plane level which is useful in
-scattering computations (OS.f,ISO.f).
+  Function: Update the atmospheric profile (P(z),T(z),H2O(z),O3(z)) in case the observer is on
+  board an aircraft.
+
+  Description: Given the altitude or pressure at aircraft level as input, the first task is to
+  compute the altitude (in case the pressure has been entered) or the pressure (in case the altitude has
+  been entered) at plane level. Then, a new atmospheric profile is created (Pp,Tp,H2Op,O3p) for which
+  the last level is located at the plane altitude. This profile is used in the gaseous absorption
+  computation (ABSTRA.f) for the path from target to sensor (upward transmission). The ozone and
+  water vapor integrated content of the "plane" atmospheric profile are also an output of this
+  subroutine. The last output is the proportion of molecules below plane level which is useful in
+  scattering computations (OS.f,ISO.f).
 */
 void Altitude::presplane(AtmosModel& atms)
 {
-	/* log linear interpolation */
-	xpp += atms.z[0];
-	if(xpp >= 100) xpp = 1000;
-
-	int i;
-	for(i = 0; atms.z[i] <= xpp; i++);
-
-	int isup = i;
-	int iinf = i-1;
-
-	float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
-	float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
-	float ps = (float)(exp((xpp - xb) / xa));
-
-	/* interpolating temperature wator vapor and ozone profile versus altitud */
-	float xalt = xpp;
-	float xtemp  = (atms.t[isup] - atms.t[iinf])/ (atms.z[isup] - atms.z[iinf]);
-	xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
-	float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
-	xwo =  xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
-	float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
-	xwh =  xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];
-
-	/* updating atmospheric profile
-	  last level: plane     , complete to 34
-	  with interpolated layers */
-	for(i = 0; i <= iinf; i++)
-	{
-		plane_sim.zpl[i] = atms.z[i];
-		plane_sim.ppl[i] = atms.p[i];
-		plane_sim.tpl[i] = atms.t[i];
-		plane_sim.whpl[i] = atms.wh[i];
-		plane_sim.wopl[i] = atms.wo[i];
-	}
+    /* log linear interpolation */
+    xpp += atms.z[0];
+    if(xpp >= 100) xpp = 1000;
+
+    int i;
+    for(i = 0; atms.z[i] <= xpp; i++);
+
+    int isup = i;
+    int iinf = i-1;
+
+    float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
+    float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
+    float ps = (float)(exp((xpp - xb) / xa));
+
+    /* interpolating temperature wator vapor and ozone profile versus altitud */
+    float xalt = xpp;
+    float xtemp  = (atms.t[isup] - atms.t[iinf])/ (atms.z[isup] - atms.z[iinf]);
+    xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
+    float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
+    xwo =  xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
+    float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
+    xwh =  xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];
+
+    /* updating atmospheric profile
+       last level: plane     , complete to 34
+       with interpolated layers */
+    for(i = 0; i <= iinf; i++)
+    {
+	plane_sim.zpl[i] = atms.z[i];
+	plane_sim.ppl[i] = atms.p[i];
+	plane_sim.tpl[i] = atms.t[i];
+	plane_sim.whpl[i] = atms.wh[i];
+	plane_sim.wopl[i] = atms.wo[i];
+    }
 
-	for(i = iinf+1; i < 34; i++)
-	{
-		plane_sim.zpl[i] = xalt;
-		plane_sim.ppl[i] = ps;
-		plane_sim.tpl[i] = xtemp;
-		plane_sim.whpl[i] = xwh;
-		plane_sim.wopl[i] = xwo;
-	}
+    for(i = iinf+1; i < 34; i++)
+    {
+	plane_sim.zpl[i] = xalt;
+	plane_sim.ppl[i] = ps;
+	plane_sim.tpl[i] = xtemp;
+	plane_sim.whpl[i] = xwh;
+	plane_sim.wopl[i] = xwo;
+    }
 
-	/* compute modified h2o and o3 integrated content
-	 compute conversion factor for rayleigh optical thickness computation
-	 ftray=rp/rt */
-	atms.uw = 0;
-	atms.uo3 = 0;
-	const float g = 98.1f;
-	const float air = 0.028964f/0.0224f;
-	const float ro3 = 0.048f/0.0224f;
-	float rt = 0;
-	float rp = 0;
-
-	float rmo3[34];
-	float rmwh[34];
-	int k;
-	for(k = 0; k < 33; k++)
-	{
-		float roair = (float)(air * 273.16 * plane_sim.ppl[k] / (1013.25 * plane_sim.tpl[k]));
-		rmwh[k] = atms.wh[k] / (roair * 1000);
-		rmo3[k] = atms.wo[k] / (roair * 1000);
-		rt += (atms.p[k+1] / atms.t[k+1] + atms.p[k] / atms.p[k]) * (atms.z[k+1] - atms.z[k]);
-		rp += (plane_sim.ppl[k+1] / plane_sim.tpl[k+1] + plane_sim.ppl[k] / plane_sim.tpl[k]) 
-			* (plane_sim.zpl[k+1] - plane_sim.zpl[k]);
-	}
+    /* compute modified h2o and o3 integrated content
+       compute conversion factor for rayleigh optical thickness computation
+       ftray=rp/rt */
+    atms.uw = 0;
+    atms.uo3 = 0;
+    const float g = 98.1f;
+    const float air = 0.028964f/0.0224f;
+    const float ro3 = 0.048f/0.0224f;
+    float rt = 0;
+    float rp = 0;
+
+    float rmo3[34];
+    float rmwh[34];
+    int k;
+    for(k = 0; k < 33; k++)
+    {
+	float roair = (float)(air * 273.16 * plane_sim.ppl[k] / (1013.25 * plane_sim.tpl[k]));
+	rmwh[k] = atms.wh[k] / (roair * 1000);
+	rmo3[k] = atms.wo[k] / (roair * 1000);
+	rt += (atms.p[k+1] / atms.t[k+1] + atms.p[k] / atms.p[k]) * (atms.z[k+1] - atms.z[k]);
+	rp += (plane_sim.ppl[k+1] / plane_sim.tpl[k+1] + plane_sim.ppl[k] / plane_sim.tpl[k]) 
+	    * (plane_sim.zpl[k+1] - plane_sim.zpl[k]);
+    }
 
-	ftray = rp / rt;
-	for(k = 1; k < 33; k++)
-	{
-		float ds = (plane_sim.ppl[k-1] - plane_sim.ppl[k]) / plane_sim.ppl[0];
-		atms.uw += (rmwh[k] + rmwh[k-1])*ds/2;
-		atms.uo3+= (rmo3[k] + rmo3[k-1])*ds/2;
-	}
+    ftray = rp / rt;
+    for(k = 1; k < 33; k++)
+    {
+	float ds = (plane_sim.ppl[k-1] - plane_sim.ppl[k]) / plane_sim.ppl[0];
+	atms.uw += (rmwh[k] + rmwh[k-1])*ds/2;
+	atms.uo3+= (rmo3[k] + rmo3[k-1])*ds/2;
+    }
 
-	atms.uw *= plane_sim.ppl[0] * 100 / g;
-	atms.uo3*= plane_sim.ppl[0] * 100 / g;
-	atms.uo3*= 1000 / ro3;
+    atms.uw *= plane_sim.ppl[0] * 100 / g;
+    atms.uo3*= plane_sim.ppl[0] * 100 / g;
+    atms.uo3*= 1000 / ro3;
 }
 
 void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
@@ -194,145 +194,145 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
     xps = original_xps;
     xpp = original_xpp;
 
-	float uwus;
-	float uo3us;
-	if(xps <= 0)
+    float uwus;
+    float uo3us;
+    if(xps <= 0)
+    {
+	xps = 0;
+	uwus = 1.424f;
+	uo3us = 0.344f;
+    }
+    else if(atms.idatm != 8) pressure(atms, atms.uw, atms.uo3);
+    else pressure(atms, uwus, uo3us);
+
+    if(xpp <= 0)
+    {
+	/* ground measurement option */
+	palt = 0;
+	pps = atms.p[0];
+	idatmp = 0;
+	taer55p = 0;
+	puw = 0;
+    }
+    else if(xpp >= 100)
+    {
+	/* satellite case of equivalent */
+	palt = 1000;
+	pps = 0;
+	taer55p = aerocon.taer55;
+	puw = 0;
+	ftray = 1;
+	idatmp = 4;
+    }
+    else
+    {
+	/* "real" plane case */
+	cin >> puw;
+	cin >> puo3;
+	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
+	if ( puw < 0 )
 	{
-		xps = 0;
-		uwus = 1.424f;
-		uo3us = 0.344f;
+	    presplane(atms);
+	    idatmp = 2;
+
+	    if (atms.idatm == 8)
+	    {
+		puwus = puw;
+		puo3us = puo3;
+		puw *= atms.uw / uwus;
+		puo3 *= atms.uo3 / uo3us;
+		idatmp = 8;
+	    }
 	}
-	else if(atms.idatm != 8) pressure(atms, atms.uw, atms.uo3);
-	else pressure(atms, uwus, uo3us);
-
-	if(xpp <= 0)
+	else
 	{
-		/* ground measurement option */
-		palt = 0;
-		pps = atms.p[0];
-		idatmp = 0;
-		taer55p = 0;
-		puw = 0;
+	    presplane(atms);
+	    idatmp = 8;
 	}
-	else if(xpp >= 100)
+
+	palt = plane_sim.zpl[33] - atms.z[0];
+	pps = plane_sim.ppl[33];
+	cin >> taer55p;
+
+	if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03))
 	{
-		/* satellite case of equivalent */
-		palt = 1000;
-		pps = 0;
-		taer55p = aerocon.taer55;
-		puw = 0;
-		ftray = 1;
-		idatmp = 4;
+	    /* a scale heigh of 2km is assumed in case no value is given for taer55p */
+	    taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 2)));
 	}
 	else
 	{
-		/* "real" plane case */
-		cin >> puw;
-		cin >> puo3;
-		cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
-		if ( puw < 0 )
-		{
-			presplane(atms);
-			idatmp = 2;
-
-			if (atms.idatm == 8)
-			{
-				puwus = puw;
-				puo3us = puo3;
-				puw *= atms.uw / uwus;
-				puo3 *= atms.uo3 / uo3us;
-				idatmp = 8;
-			}
-		}
-		else
-		{
-			presplane(atms);
-			idatmp = 8;
-		}
-
-		palt = plane_sim.zpl[33] - atms.z[0];
-		pps = plane_sim.ppl[33];
-		cin >> taer55p;
-
-		if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03))
-		{
-			/* a scale heigh of 2km is assumed in case no value is given for taer55p */
-			taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 2)));
-		}
-		else
-		{
-			/* compute effective scale heigh */
-			double sham = exp(-palt / 4);
-			double sha = 1 - (taer55p / aerocon.taer55);
-
-			if( sha >= sham) taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
-			else {
-				sha = -palt/log(sha);
-				taer55p = (float)(aerocon.taer55 * (1 - exp(-palt/sha)));
-			}
-		}
+	    /* compute effective scale heigh */
+	    double sham = exp(-palt / 4);
+	    double sha = 1 - (taer55p / aerocon.taer55);
+
+	    if( sha >= sham) taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
+	    else {
+		sha = -palt/log(sha);
+		taer55p = (float)(aerocon.taer55 * (1 - exp(-palt/sha)));
+	    }
 	}
+    }
 }
 
 void Altitude::parse()
 {
-	cin >> original_xps;
-	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
+    cin >> original_xps;
+    cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
     original_xps = -original_xps;
     
-	cin >> original_xpp;
-	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
+    cin >> original_xpp;
+    cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
     original_xpp = -original_xpp;
 }
 
 /* --- plane simulation output if selected ---- */
 void Altitude::print()
 {
-	if(palt < 1000)
-	{
-		Output::Ln();
-		Output::WriteLn(22," plane simulation description ");
-		Output::WriteLn(22," ---------------------------- ");
+    if(palt < 1000)
+    {
+	Output::Ln();
+	Output::WriteLn(22," plane simulation description ");
+	Output::WriteLn(22," ---------------------------- ");
 		
-		ostringstream s1;
-		s1.setf(ios::fixed, ios::floatfield);
-		s1.precision(2);
-		s1 << " plane  pressure          [mb] " << setw(9) << pps << ends;
-		Output::WriteLn(10,s1.str());
-
-		ostringstream s2;
-		s2.setf(ios::fixed, ios::floatfield);
-		s2.precision(3);
-		s2 << " plane  altitude absolute [km] " << setw(9) << plane_sim.zpl[33] << ends;
-		Output::WriteLn(10,s2.str());
+	ostringstream s1;
+	s1.setf(ios::fixed, ios::floatfield);
+	s1.precision(2);
+	s1 << " plane  pressure          [mb] " << setw(9) << pps << ends;
+	Output::WriteLn(10,s1.str());
+
+	ostringstream s2;
+	s2.setf(ios::fixed, ios::floatfield);
+	s2.precision(3);
+	s2 << " plane  altitude absolute [km] " << setw(9) << plane_sim.zpl[33] << ends;
+	Output::WriteLn(10,s2.str());
 
 		
-		Output::WriteLn(15," atmosphere under plane description: ");
-
-		ostringstream s3;
-		s3.setf(ios::fixed, ios::floatfield);
-		s3.precision(3);
-		s3 << " ozone content            " << setw(9) << puo3 << ends;
-		Output::WriteLn(15,s3.str());
-
-
-		ostringstream s4;
-		s4.setf(ios::fixed, ios::floatfield);
-		s4.precision(3);
-		s4 << " h2o   content            " << setw(9) << puw << ends;
-		Output::WriteLn(15,s4.str());
-
-		ostringstream s5;
-		s5.setf(ios::fixed, ios::floatfield);
-		s5.precision(3);
-		s5 << "aerosol opt. thick. 550nm " << setw(9) << taer55p << ends;
-		Output::WriteLn(15,s5.str());
-	}
+	Output::WriteLn(15," atmosphere under plane description: ");
+
+	ostringstream s3;
+	s3.setf(ios::fixed, ios::floatfield);
+	s3.precision(3);
+	s3 << " ozone content            " << setw(9) << puo3 << ends;
+	Output::WriteLn(15,s3.str());
+
+
+	ostringstream s4;
+	s4.setf(ios::fixed, ios::floatfield);
+	s4.precision(3);
+	s4 << " h2o   content            " << setw(9) << puw << ends;
+	Output::WriteLn(15,s4.str());
+
+	ostringstream s5;
+	s5.setf(ios::fixed, ios::floatfield);
+	s5.precision(3);
+	s5 << "aerosol opt. thick. 550nm " << setw(9) << taer55p << ends;
+	Output::WriteLn(15,s5.str());
+    }
 }
 
 Altitude Altitude::Parse()
 {
-	Altitude alt;
-	alt.parse();
-	return alt;
+    Altitude alt;
+    alt.parse();
+    return alt;
 }

imagery/i.atcorr/AtmosModel.cpp

@@ -3,315 +3,315 @@
 
 void AtmosModel::tropic()
 {
-	static const float z1[34] =
+    static const float z1[34] =
 	{ 
-		0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f, 
-		12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 
-		22.f, 23.f, 24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
+	    0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f, 
+	    12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 
+	    22.f, 23.f, 24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
 	};
 	
-	static const float p1[34] =
-    { 
-		1013.f, 904.f, 805.f, 715.f, 633.f, 559.f, 492.f, 432.f, 378.f, 
-		329.f, 286.f, 247.f, 213.f, 182.f, 156.f, 132.f, 111.f, 93.7f,
-		78.9f, 66.6f, 56.5f, 48.f, 40.9f, 35.f, 30.f, 25.7f, 12.2f, 6.f, 
-		3.05f, 1.59f, .854f, .0579f, 3e-4f, 0.f
+    static const float p1[34] =
+	{ 
+	    1013.f, 904.f, 805.f, 715.f, 633.f, 559.f, 492.f, 432.f, 378.f, 
+	    329.f, 286.f, 247.f, 213.f, 182.f, 156.f, 132.f, 111.f, 93.7f,
+	    78.9f, 66.6f, 56.5f, 48.f, 40.9f, 35.f, 30.f, 25.7f, 12.2f, 6.f, 
+	    3.05f, 1.59f, .854f, .0579f, 3e-4f, 0.f
 	};
 
-	static const float t1[34] =
-    { 
-		300.f, 294.f, 288.f, 284.f, 277.f, 270.f, 264.f, 257.f, 250.f, 
-		244.f, 237.f, 230.f, 224.f, 217.f, 210.f, 204.f, 197.f, 195.f,
-		199.f, 203.f, 207.f, 211.f, 215.f, 217.f, 219.f, 221.f, 232.f, 
-		243.f, 254.f, 265.f, 270.f, 219.f, 210.f, 210.f
+    static const float t1[34] =
+	{ 
+	    300.f, 294.f, 288.f, 284.f, 277.f, 270.f, 264.f, 257.f, 250.f, 
+	    244.f, 237.f, 230.f, 224.f, 217.f, 210.f, 204.f, 197.f, 195.f,
+	    199.f, 203.f, 207.f, 211.f, 215.f, 217.f, 219.f, 221.f, 232.f, 
+	    243.f, 254.f, 265.f, 270.f, 219.f, 210.f, 210.f
 	};
 
-	static const float wh1[34] =
-    { 
-		19.f, 13.f, 9.3f, 4.7f, 2.2f, 1.5f, .85f, .47f, .25f, .12f, .05f, 
+    static const float wh1[34] =
+	{ 
+	    19.f, 13.f, 9.3f, 4.7f, 2.2f, 1.5f, .85f, .47f, .25f, .12f, .05f, 
 
 
-		.017f, .006f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f, 5e-4f,
-		4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f, 
-		3.6e-4f, 1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
+	    .017f, .006f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f, 5e-4f,
+	    4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f, 
+	    3.6e-4f, 1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
 	};
 
-	static const float wo1[34] =
-    { 
-		5.6e-5f, 5.6e-5f, 5.4e-5f, 5.1e-5f, 4.7e-5f, 4.5e-5f,
-		4.3e-5f, 4.1e-5f, 3.9e-5f, 3.9e-5f, 3.9e-5f, 4.1e-5f, 4.3e-5f, 4.5e-5f,
+    static const float wo1[34] =
+	{ 
+	    5.6e-5f, 5.6e-5f, 5.4e-5f, 5.1e-5f, 4.7e-5f, 4.5e-5f,
+	    4.3e-5f, 4.1e-5f, 3.9e-5f, 3.9e-5f, 3.9e-5f, 4.1e-5f, 4.3e-5f, 4.5e-5f,
 	    4.5e-5f, 4.7e-5f, 4.7e-5f, 6.9e-5f, 9e-5f, 1.4e-4f, 1.9e-4f, 2.4e-4f,
-		2.8e-4f, 3.2e-4f, 3.4e-4f, 3.4e-4f, 2.4e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f,
-		4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
+	    2.8e-4f, 3.2e-4f, 3.4e-4f, 3.4e-4f, 2.4e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f,
+	    4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
 	};
 
-	/* model: tropical mc clatchey */
-	for (int i = 0; i < 34; i++)
+    /* model: tropical mc clatchey */
+    for (int i = 0; i < 34; i++)
     {
-      z[i] = z1[i];
-      p[i] = p1[i];
-      t[i] = t1[i];
-      wh[i] = wh1[i];
-      wo[i] = wo1[i];
+	z[i] = z1[i];
+	p[i] = p1[i];
+	t[i] = t1[i];
+	wh[i] = wh1[i];
+	wo[i] = wo1[i];
     }
 }
 
 void AtmosModel::midsum()
 {
-	static const float z1[34] =
-    { 
-		0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
-		12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
-		24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
+    static const float z1[34] =
+	{ 
+	    0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
+	    12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
+	    24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
 	};
 
-	static const float p1[34] =
+    static const float p1[34] =
 	{ 
-		1013.f, 902.f, 802.f, 710.f, 628.f, 554.f, 487.f, 426.f,
-		372.f, 324.f, 281.f, 243.f, 209.f, 179.f, 153.f, 130.f, 111.f, 95.f,
-		81.2f, 69.5f, 59.5f, 51.f, 43.7f, 37.6f, 32.2f, 27.7f, 13.2f, 6.52f, 
-		3.33f, 1.76f, .951f, .0671f, 3e-4f, 0.f
+	    1013.f, 902.f, 802.f, 710.f, 628.f, 554.f, 487.f, 426.f,
+	    372.f, 324.f, 281.f, 243.f, 209.f, 179.f, 153.f, 130.f, 111.f, 95.f,
+	    81.2f, 69.5f, 59.5f, 51.f, 43.7f, 37.6f, 32.2f, 27.7f, 13.2f, 6.52f, 
+	    3.33f, 1.76f, .951f, .0671f, 3e-4f, 0.f
 	};
 
-	static const float t1[34] =
-    { 
-		294.f, 290.f, 285.f, 279.f, 273.f, 267.f, 261.f, 255.f,
+    static const float t1[34] =
+	{ 
+	    294.f, 290.f, 285.f, 279.f, 273.f, 267.f, 261.f, 255.f,
 	    248.f, 242.f, 235.f, 229.f, 222.f, 216.f, 216.f, 216.f, 216.f, 216.f,
-		216.f, 217.f, 218.f, 219.f, 220.f, 222.f, 223.f, 224.f, 234.f, 245.f, 258.f,
-		270.f, 276.f, 218.f, 210.f, 210.f
+	    216.f, 217.f, 218.f, 219.f, 220.f, 222.f, 223.f, 224.f, 234.f, 245.f, 258.f,
+	    270.f, 276.f, 218.f, 210.f, 210.f
 	};
 
-	static const float wh1[34] =
-    { 
-		14.f, 9.3f, 5.9f, 3.3f, 1.9f, 1.f, .61f, .37f, .21f, .12f,
-		.064f, .022f, .006f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f, 5e-4f,
-		4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f, 3.6e-4f,
-		1.1e-4f, 4.3e-5f, 1.9e-5f, 1.3e-6f, 1.4e-7f, 1e-9f, 0.f
+    static const float wh1[34] =
+	{ 
+	    14.f, 9.3f, 5.9f, 3.3f, 1.9f, 1.f, .61f, .37f, .21f, .12f,
+	    .064f, .022f, .006f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f, 5e-4f,
+	    4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f, 3.6e-4f,
+	    1.1e-4f, 4.3e-5f, 1.9e-5f, 1.3e-6f, 1.4e-7f, 1e-9f, 0.f
 	};
 
-	static const float wo1[34] =
-    { 
-		6e-5f, 6e-5f, 6e-5f, 6.2e-5f, 6.4e-5f, 6.6e-5f, 6.9e-5f,
+    static const float wo1[34] =
+	{ 
+	    6e-5f, 6e-5f, 6e-5f, 6.2e-5f, 6.4e-5f, 6.6e-5f, 6.9e-5f,
 	    7.5e-5f, 7.9e-5f, 8.6e-5f, 9e-5f, 1.1e-4f, 1.2e-4f, 1.5e-4f, 1.8e-4f,
-		1.9e-4f, 2.1e-4f, 2.4e-4f, 2.8e-4f, 3.2e-4f, 3.4e-4f, 3.6e-4f, 3.6e-4f,
-		3.4e-4f, 3.2e-4f, 3e-4f, 2e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f, 4.3e-6f,
-		8.6e-8f, 4.3e-11f, 0.f
+	    1.9e-4f, 2.1e-4f, 2.4e-4f, 2.8e-4f, 3.2e-4f, 3.4e-4f, 3.6e-4f, 3.6e-4f,
+	    3.4e-4f, 3.2e-4f, 3e-4f, 2e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f, 4.3e-6f,
+	    8.6e-8f, 4.3e-11f, 0.f
 	};
 
-	/* model: midlatitude summer mc clatchey */
-	for (int i = 0; i < 34; i++)
+    /* model: midlatitude summer mc clatchey */
+    for (int i = 0; i < 34; i++)
     {
-      z[i] = z1[i];
-      p[i] = p1[i];
-      t[i] = t1[i];
-      wh[i] = wh1[i];
-      wo[i] = wo1[i];
+	z[i] = z1[i];
+	p[i] = p1[i];
+	t[i] = t1[i];
+	wh[i] = wh1[i];
+	wo[i] = wo1[i];
     }
 }
 
 void AtmosModel::midwin()
 {
-	static const float z1[34] =
+    static const float z1[34] =
 	{ 
-		0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
-		12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
-		24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
+	    0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
+	    12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
+	    24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
 	};
 
-	static const float p1[34] =
-    { 
-		1018.f, 897.3f, 789.7f, 693.8f, 608.1f, 531.3f, 462.7f,
+    static const float p1[34] =
+	{ 
+	    1018.f, 897.3f, 789.7f, 693.8f, 608.1f, 531.3f, 462.7f,
 	    401.6f, 347.3f, 299.2f, 256.8f, 219.9f, 188.2f, 161.f, 137.8f, 117.8f,
-		100.7f, 86.1f, 73.5f, 62.8f, 53.7f, 45.8f, 39.1f, 33.4f, 28.6f, 24.3f,
+	    100.7f, 86.1f, 73.5f, 62.8f, 53.7f, 45.8f, 39.1f, 33.4f, 28.6f, 24.3f,
 	    11.1f, 5.18f, 2.53f, 1.29f, .682f, .0467f, 3e-4f, 0.f
 	};
 
-	static const float t1[34] =
-    { 
-		272.2f, 268.7f, 265.2f, 261.7f, 255.7f, 249.7f, 243.7f,
+    static const float t1[34] =
+	{ 
+	    272.2f, 268.7f, 265.2f, 261.7f, 255.7f, 249.7f, 243.7f,
 	    237.7f, 231.7f, 225.7f, 219.7f, 219.2f, 218.7f, 218.2f, 217.7f, 217.2f,
-		216.7f, 216.2f, 215.7f, 215.2f, 215.2f, 215.2f, 215.2f, 215.2f, 215.2f,
+	    216.7f, 216.2f, 215.7f, 215.2f, 215.2f, 215.2f, 215.2f, 215.2f, 215.2f,
 	    215.2f, 217.4f, 227.8f, 243.2f, 258.5f, 265.7f, 230.7f, 210.2f, 210.f
 	};	
 
-	static const float wh1[34] =
-    { 
-		3.5f, 2.5f, 1.8f, 1.2f, .66f, .38f, .21f, .085f, .035f,
+    static const float wh1[34] =
+	{ 
+	    3.5f, 2.5f, 1.8f, 1.2f, .66f, .38f, .21f, .085f, .035f,
 	    .016f, .0075f, .0069f, .006f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f,
-		5e-4f, 4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f,
+	    5e-4f, 4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f,
 	    3.6e-4f, 1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
 	};
 
-	static const float wo1[34] = 
+    static const float wo1[34] = 
 	{ 
-		6e-5f, 5.4e-5f, 4.9e-5f, 4.9e-5f, 4.9e-5f, 5.8e-5f,
+	    6e-5f, 5.4e-5f, 4.9e-5f, 4.9e-5f, 4.9e-5f, 5.8e-5f,
 	    6.4e-5f, 7.7e-5f, 9e-5f, 1.2e-4f, 1.6e-4f, 2.1e-4f, 2.6e-4f, 3e-4f,
-		3.2e-4f, 3.4e-4f, 3.6e-4f, 3.9e-4f, 4.1e-4f, 4.3e-4f, 4.5e-4f, 4.3e-4f,
+	    3.2e-4f, 3.4e-4f, 3.6e-4f, 3.9e-4f, 4.1e-4f, 4.3e-4f, 4.5e-4f, 4.3e-4f,
 	    4.3e-4f, 3.9e-4f, 3.6e-4f, 3.4e-4f, 1.9e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f,
-		4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
+	    4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
 	};
 
-	/* model: midlatitude winter mc clatchey */
-	for (int i = 0; i < 34; i++)
+    /* model: midlatitude winter mc clatchey */
+    for (int i = 0; i < 34; i++)
     {
-      z[i] = z1[i];
-      p[i] = p1[i];
-      t[i] = t1[i];
-      wh[i] = wh1[i];
-      wo[i] = wo1[i];
+	z[i] = z1[i];
+	p[i] = p1[i];
+	t[i] = t1[i];
+	wh[i] = wh1[i];
+	wo[i] = wo1[i];
     }
 }
 
 void AtmosModel::subsum()
 {
-	static const float z1[34] =
-    { 
-		0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
+    static const float z1[34] =
+	{ 
+	    0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
 	    12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
-		24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
+	    24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
 	};
 
-	static const float p1[34] =
-    { 
-		1010.f, 896.f, 792.9f, 700.f, 616.f, 541.f, 473.f, 413.f,
+    static const float p1[34] =
+	{ 
+	    1010.f, 896.f, 792.9f, 700.f, 616.f, 541.f, 473.f, 413.f,
 	    359.f, 310.7f, 267.7f, 230.f, 197.7f, 170.f, 146.f, 125.f, 108.f, 92.8f,
-		79.8f, 68.6f, 58.9f, 50.7f, 43.6f, 37.5f, 32.27f, 27.8f, 13.4f, 6.61f,
+	    79.8f, 68.6f, 58.9f, 50.7f, 43.6f, 37.5f, 32.27f, 27.8f, 13.4f, 6.61f,
 	    3.4f, 1.81f, .987f, .0707f, 3e-4f, 0.f
 	};
 
-	static const float t1[34] =
-    { 
-		287.f, 282.f, 276.f, 271.f, 266.f, 260.f, 253.f, 246.f,
-		239.f, 232.f, 225.f, 225.f, 225.f, 225.f, 225.f, 225.f, 225.f, 225.f,
+    static const float t1[34] =
+	{ 
+	    287.f, 282.f, 276.f, 271.f, 266.f, 260.f, 253.f, 246.f,
+	    239.f, 232.f, 225.f, 225.f, 225.f, 225.f, 225.f, 225.f, 225.f, 225.f,
 	    225.f, 225.f, 225.f, 225.f, 225.f, 225.f, 226.f, 228.f, 235.f, 247.f, 262.f,
-		274.f, 277.f, 216.f, 210.f, 210.f
+	    274.f, 277.f, 216.f, 210.f, 210.f
 	};
 
-	static const float wh1[34] =
-    { 
-		9.1f, 6.f, 4.2f, 2.7f, 1.7f, 1.f, .54f, .29f, .13f, .042f,
+    static const float wh1[34] =
+	{ 
+	    9.1f, 6.f, 4.2f, 2.7f, 1.7f, 1.f, .54f, .29f, .13f, .042f,
 	    .015f, .0094f, .006f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f, 5e-4f,
-		4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f, 3.6e-4f,
-		1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
+	    4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f, 3.6e-4f,
+	    1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
 	};
 
-	static const float wo1[34] = 
+    static const float wo1[34] = 
 	{ 
-		4.9e-5f, 5.4e-5f, 5.6e-5f, 5.8e-5f, 6e-5f, 6.4e-5f,
+	    4.9e-5f, 5.4e-5f, 5.6e-5f, 5.8e-5f, 6e-5f, 6.4e-5f,
 	    7.1e-5f, 7.5e-5f, 7.9e-5f, 1.1e-4f, 1.3e-4f, 1.8e-4f, 2.1e-4f, 2.6e-4f,
-		2.8e-4f, 3.2e-4f, 3.4e-4f, 3.9e-4f, 4.1e-4f, 4.1e-4f, 3.9e-4f, 3.6e-4f,
+	    2.8e-4f, 3.2e-4f, 3.4e-4f, 3.9e-4f, 4.1e-4f, 4.1e-4f, 3.9e-4f, 3.6e-4f,
 	    3.2e-4f, 3e-4f, 2.8e-4f, 2.6e-4f, 1.4e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f,
-		4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
+	    4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
 	};
 
-	/* model: subarctique summer mc clatchey */
-	for (int i = 0; i < 34; i++)
+    /* model: subarctique summer mc clatchey */
+    for (int i = 0; i < 34; i++)
     {
-      z[i] = z1[i];
-      p[i] = p1[i];
-      t[i] = t1[i];
-      wh[i] = wh1[i];
-      wo[i] = wo1[i];
+	z[i] = z1[i];
+	p[i] = p1[i];
+	t[i] = t1[i];
+	wh[i] = wh1[i];
+	wo[i] = wo1[i];
     }
 }
 
 void AtmosModel::subwin()
 {
-	static const float z1[34] =
-    { 
-		0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
+    static const float z1[34] =
+	{ 
+	    0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
 	    12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
-		24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
+	    24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
 	};
 
-	static const float p1[34] =
-    { 
-		1013.f, 887.8f, 777.5f, 679.8f, 593.2f, 515.8f, 446.7f,
+    static const float p1[34] =
+	{ 
+	    1013.f, 887.8f, 777.5f, 679.8f, 593.2f, 515.8f, 446.7f,
 	    385.3f, 330.8f, 282.9f, 241.8f, 206.7f, 176.6f, 151.f, 129.1f, 110.3f,
-		94.31f, 80.58f, 68.82f, 58.75f, 50.14f, 42.77f, 36.47f, 31.09f, 26.49f,
-		22.56f, 10.2f, 4.701f, 2.243f, 1.113f, .5719f, .04016f, 3e-4f, 0.f
+	    94.31f, 80.58f, 68.82f, 58.75f, 50.14f, 42.77f, 36.47f, 31.09f, 26.49f,
+	    22.56f, 10.2f, 4.701f, 2.243f, 1.113f, .5719f, .04016f, 3e-4f, 0.f
 	};
 
-	static const float t1[34] =
-    { 
-		257.1f, 259.1f, 255.9f, 252.7f, 247.7f, 240.9f, 234.1f,
+    static const float t1[34] =
+	{ 
+	    257.1f, 259.1f, 255.9f, 252.7f, 247.7f, 240.9f, 234.1f,
 	    227.3f, 220.6f, 217.2f, 217.2f, 217.2f, 217.2f, 217.2f, 217.2f, 217.2f,
-		216.6f, 216.f, 215.4f, 214.8f, 214.1f, 213.6f, 213.f, 212.4f, 211.8f,
-		211.2f, 216.f, 222.2f, 234.7f, 247.f, 259.3f, 245.7f, 210.f, 210.f
+	    216.6f, 216.f, 215.4f, 214.8f, 214.1f, 213.6f, 213.f, 212.4f, 211.8f,
+	    211.2f, 216.f, 222.2f, 234.7f, 247.f, 259.3f, 245.7f, 210.f, 210.f
 	};
 
-	static const float wh1[34] =
-    { 
-		1.2f, 1.2f, .94f, .68f, .41f, .2f, .098f, .054f, .011f,
-		.0084f, .0055f, .0038f, .0026f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f,
-		5e-4f, 4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f,
-		3.6e-4f, 1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
+    static const float wh1[34] =
+	{ 
+	    1.2f, 1.2f, .94f, .68f, .41f, .2f, .098f, .054f, .011f,
+	    .0084f, .0055f, .0038f, .0026f, .0018f, .001f, 7.6e-4f, 6.4e-4f, 5.6e-4f,
+	    5e-4f, 4.9e-4f, 4.5e-4f, 5.1e-4f, 5.1e-4f, 5.4e-4f, 6e-4f, 6.7e-4f,
+	    3.6e-4f, 1.1e-4f, 4.3e-5f, 1.9e-5f, 6.3e-6f, 1.4e-7f, 1e-9f, 0.f
 	};
 
-	static const float wo1[34] =
-    { 
-		4.1e-5f, 4.1e-5f, 4.1e-5f, 4.3e-5f, 4.5e-5f, 4.7e-5f,
-		4.9e-5f, 7.1e-5f, 9e-5f, 1.6e-4f, 2.4e-4f, 3.2e-4f, 4.3e-4f, 4.7e-4f,
-		4.9e-4f, 5.6e-4f, 6.2e-4f, 6.2e-4f, 6.2e-4f, 6e-4f, 5.6e-4f, 5.1e-4f,
-		4.7e-4f, 4.3e-4f, 3.6e-4f, 3.2e-4f, 1.5e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f,
-		4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
+    static const float wo1[34] =
+	{ 
+	    4.1e-5f, 4.1e-5f, 4.1e-5f, 4.3e-5f, 4.5e-5f, 4.7e-5f,
+	    4.9e-5f, 7.1e-5f, 9e-5f, 1.6e-4f, 2.4e-4f, 3.2e-4f, 4.3e-4f, 4.7e-4f,
+	    4.9e-4f, 5.6e-4f, 6.2e-4f, 6.2e-4f, 6.2e-4f, 6e-4f, 5.6e-4f, 5.1e-4f,
+	    4.7e-4f, 4.3e-4f, 3.6e-4f, 3.2e-4f, 1.5e-4f, 9.2e-5f, 4.1e-5f, 1.3e-5f,
+	    4.3e-6f, 8.6e-8f, 4.3e-11f, 0.f
 	};
 
-	/* model: subarctique winter mc clatchey */
-	for (int i = 0; i < 34; i++)
+    /* model: subarctique winter mc clatchey */
+    for (int i = 0; i < 34; i++)
     {
-      z[i] = z1[i];
-      p[i] = p1[i];
-      t[i] = t1[i];
-      wh[i] = wh1[i];
-      wo[i] = wo1[i];
+	z[i] = z1[i];
+	p[i] = p1[i];
+	t[i] = t1[i];
+	wh[i] = wh1[i];
+	wo[i] = wo1[i];
     }
 }
 
 void AtmosModel::us62()
 {
-	static const float z1[34] =
-    { 
-		0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
+    static const float z1[34] =
+	{ 
+	    0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 11.f,
 	    12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f, 20.f, 21.f, 22.f, 23.f,
-		24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
+	    24.f, 25.f, 30.f, 35.f, 40.f, 45.f, 50.f, 70.f, 100.f, 99999.f
 	};
 
-	static const float p1[34] =
-    { 
-		1013.f, 898.6f, 795.f, 701.2f, 616.6f, 540.5f, 472.2f,
-		411.1f, 356.5f, 308.f, 265.f, 227.f, 194.f, 165.8f, 141.7f, 121.1f,
-		103.5f, 88.5f, 75.65f, 64.67f, 55.29f, 47.29f, 40.47f, 34.67f, 29.72f, 25.49f,
-		11.97f, 5.746f, 2.871f, 1.491f, .7978f, .0552f, 3.008e-4f, 0.f
+    static const float p1[34] =
+	{ 
+	    1013.f, 898.6f, 795.f, 701.2f, 616.6f, 540.5f, 472.2f,
+	    411.1f, 356.5f, 308.f, 265.f, 227.f, 194.f, 165.8f, 141.7f, 121.1f,
+	    103.5f, 88.5f, 75.65f, 64.67f, 55.29f, 47.29f, 40.47f, 34.67f, 29.72f, 25.49f,
+	    11.97f, 5.746f, 2.871f, 1.491f, .7978f, .0552f, 3.008e-4f, 0.f
 	};
 
-	static const float t1[34] =
-    { 
-		288.1f, 281.6f, 275.1f, 268.7f, 262.2f, 255.7f, 249.2f,
+    static const float t1[34] =
+	{ 
+	    288.1f, 281.6f, 275.1f, 268.7f, 262.2f, 255.7f, 249.2f,
 	    242.7f, 236.2f, 229.7f, 223.2f, 216.8f, 216.6f, 216.6f, 216.6f, 216.6f,
-		216.6f, 216.6f, 216.6f, 216.6f, 216.6f, 217.6f, 218.6f, 219.6f, 220.6f,
-		221.6f, 226.5f, 236.5f, 253.4f, 264.2f, 270.6f, 219.7f, 210.f, 210.f
+	    216.6f, 216.6f, 216.6f, 216.6f, 216.6f, 217.6f, 218.6f, 219.6f, 220.6f,
+	    221.6f, 226.5f, 236.5f, 253.4f, 264.2f, 270.6f, 219.7f, 210.f, 210.f
 	};
 
-	static const float wh1[34] =
-    { 
-		5.9f, 4.2f, 2.9f, 1.8f, 1.1f, .64f, .38f, .21f, .12f,
+    static const float wh1[34] =
+	{ 
+	    5.9f, 4.2f, 2.9f, 1.8f, 1.1f, .64f, .38f, .21f, .12f,
 	    .046f, .018f, .0082f, .0037f, .0018f, 8.4e-4f, 7.2e-4f, 6.1e-4f, 5.2e-4f,
-		4.4e-4f, 4.4e-4f, 4.4e-4f, 4.8e-4f, 5.2e-4f, 5.7e-4f, 6.1e-4f, 6.6e-4f,
+	    4.4e-4f, 4.4e-4f, 4.4e-4f, 4.8e-4f, 5.2e-4f, 5.7e-4f, 6.1e-4f, 6.6e-4f,
 	    3.8e-4f, 1.6e-4f, 6.7e-5f, 3.2e-5f, 1.2e-5f, 1.5e-7f, 1e-9f, 0.f
 	};
 	
-	static const float wo1[34] = 
+    static const float wo1[34] = 
 	{ 
-		5.4e-5f, 5.4e-5f, 5.4e-5f, 5e-5f, 4.6e-5f, 4.6e-5f,
-		4.5e-5f, 4.9e-5f, 5.2e-5f, 7.1e-5f, 9e-5f, 1.3e-4f, 1.6e-4f, 1.7e-4f,
-		1.9e-4f, 2.1e-4f, 2.4e-4f, 2.8e-4f, 3.2e-4f, 3.5e-4f, 3.8e-4f, 3.8e-4f,
-		3.9e-4f, 3.8e-4f, 3.6e-4f, 3.4e-4f, 2e-4f, 1.1e-4f, 4.9e-5f, 1.7e-5f,
-		4e-6f, 8.6e-8f, 4.3e-11f, 0.f
+	    5.4e-5f, 5.4e-5f, 5.4e-5f, 5e-5f, 4.6e-5f, 4.6e-5f,
+	    4.5e-5f, 4.9e-5f, 5.2e-5f, 7.1e-5f, 9e-5f, 1.3e-4f, 1.6e-4f, 1.7e-4f,
+	    1.9e-4f, 2.1e-4f, 2.4e-4f, 2.8e-4f, 3.2e-4f, 3.5e-4f, 3.8e-4f, 3.8e-4f,
+	    3.9e-4f, 3.8e-4f, 3.6e-4f, 3.4e-4f, 2e-4f, 1.1e-4f, 4.9e-5f, 1.7e-5f,
+	    4e-6f, 8.6e-8f, 4.3e-11f, 0.f
 	};
 
-	/* model: us standard 62 mc clatchey */
-	for (int i = 0; i < 34; i++)
+    /* model: us standard 62 mc clatchey */
+    for (int i = 0; i < 34; i++)
     {
         z[i] = z1[i];
         p[i] = p1[i];
@@ -324,142 +324,142 @@ void AtmosModel::us62()
 
 void AtmosModel::parse()
 {
-	cin >> idatm;
-	cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
+    cin >> idatm;
+    cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
 
-	uw = 0.;
-	uo3 = 0.;
+    uw = 0.;
+    uo3 = 0.;
 
-	switch(idatm)
+    switch(idatm)
+    {
+    case 0: us62();	    break;
+    case 1: tropic();	break;
+    case 2: midsum();	break;
+    case 3: midwin();	break; 
+    case 4: subsum();	break;
+    case 5: subwin();	break;
+    case 6: us62();	    break;
+    case 7: 
+    {
+	/* read input */
+	for(int i = 0; i < 34; i++)
 	{
-	case 0: us62();	    break;
-	case 1: tropic();	break;
-	case 2: midsum();	break;
-	case 3: midwin();	break; 
-	case 4: subsum();	break;
-	case 5: subwin();	break;
-	case 6: us62();	    break;
-	case 7: 
-		{
-			/* read input */
-			for(int i = 0; i < 34; i++)
-			{
-				cin >> z[i];
-				cin >> p[i];
-				cin >> t[i];
-				cin >> wh[i];
-				cin >> wo[i];
-				cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
-			}
-			break;
-		}
-	case 8: 
-		{
-			cin >> uw;
-			cin >> uo3;
-			cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
-			us62();
-            break;
-		}
-	default: fprintf(stderr, "Unknown atmospheric model!\n");
+	    cin >> z[i];
+	    cin >> p[i];
+	    cin >> t[i];
+	    cin >> wh[i];
+	    cin >> wo[i];
+	    cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
 	}
+	break;
+    }
+    case 8: 
+    {
+	cin >> uw;
+	cin >> uo3;
+	cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
+	us62();
+	break;
+    }
+    default: fprintf(stderr, "Unknown atmospheric model!\n");
+    }
 }
 
 /* --- atmospheric model ---- */
 void AtmosModel::print()
 {	
-	static const string head(" atmospheric model description  ");
-	static const string line(" -----------------------------  ");
-	Output::Begin(); Output::Repeat(22,' '); Output::Print(head); Output::End();
-	Output::Begin(); Output::Repeat(22,' '); Output::Print(line); Output::End();
+    static const string head(" atmospheric model description  ");
+    static const string line(" -----------------------------  ");
+    Output::Begin(); Output::Repeat(22,' '); Output::Print(head); Output::End();
+    Output::Begin(); Output::Repeat(22,' '); Output::Print(line); Output::End();
 
-	if(idatm < 7) 
-	{
-		static const string atmid[7] = {
-			string("no absorption computed                             "),
-			string("tropical            (uh2o=4.12g/cm2,uo3=.247cm-atm)"),
-			string("midlatitude summer  (uh2o=2.93g/cm2,uo3=.319cm-atm)"),
-			string("midlatitude winter  (uh2o=.853g/cm2,uo3=.395cm-atm)"),
-			string("subarctic  summer   (uh2o=2.10g/cm2,uo3=.480cm-atm)"),
-			string("subarctic  winter   (uh2o=.419g/cm2,uo3=.480cm-atm)"),
-			string("us  standard 1962   (uh2o=1.42g/cm2,uo3=.344cm-atm)")
-		};
-
-		Output::Begin(); 
-		Output::Repeat(10,' ');
-		Output::Print(" atmospheric model identity : ");
-		Output::End();
-
-		Output::Begin(); 
-		Output::Repeat(15,' ');
-		Output::Print(atmid[idatm]);
-		Output::End();
-	}
-	else if(idatm == 7)
-	{
-		Output::Begin();
-		Output::Print(" atmospheric model identity : ");
-		Output::End();
-
-		Output::Begin();
-		Output::Repeat(12, ' ');
-		Output::Print(" user defined atmospheric model  ");
-		Output::End();
-
-		Output::Begin();
-		Output::Repeat(12, ' ');
-		Output::Print("*altitude  *pressure  *temp.     *h2o dens. *o3 dens.  ");
-		Output::End();
-
-		for(int i = 0; i < 34; i++)
-		{
-			Output::Begin();
-			Output::Repeat(12, ' ');
-			ostringstream s;
-			s.setf(ios::fixed, ios::floatfield);
-			s << setprecision(4);
-			s << setw(9) << z[i] << "  ";
-			s << setw(9) << p[i] << "  ";
-			s << setw(9) << t[i] << "  ";
-			s << setw(9) << wh[i] << "  ";
-			s << setw(9) << wo[i] << "  ";
-			s << ends;
-			Output::Print(s.str());
-			Output::End();
-		}
-	}
-	else 
+    if(idatm < 7) 
+    {
+	static const string atmid[7] = {
+	    string("no absorption computed                             "),
+	    string("tropical            (uh2o=4.12g/cm2,uo3=.247cm-atm)"),
+	    string("midlatitude summer  (uh2o=2.93g/cm2,uo3=.319cm-atm)"),
+	    string("midlatitude winter  (uh2o=.853g/cm2,uo3=.395cm-atm)"),
+	    string("subarctic  summer   (uh2o=2.10g/cm2,uo3=.480cm-atm)"),
+	    string("subarctic  winter   (uh2o=.419g/cm2,uo3=.480cm-atm)"),
+	    string("us  standard 1962   (uh2o=1.42g/cm2,uo3=.344cm-atm)")
+	};
+
+	Output::Begin(); 
+	Output::Repeat(10,' ');
+	Output::Print(" atmospheric model identity : ");
+	Output::End();
+
+	Output::Begin(); 
+	Output::Repeat(15,' ');
+	Output::Print(atmid[idatm]);
+	Output::End();
+    }
+    else if(idatm == 7)
+    {
+	Output::Begin();
+	Output::Print(" atmospheric model identity : ");
+	Output::End();
+
+	Output::Begin();
+	Output::Repeat(12, ' ');
+	Output::Print(" user defined atmospheric model  ");
+	Output::End();
+
+	Output::Begin();
+	Output::Repeat(12, ' ');
+	Output::Print("*altitude  *pressure  *temp.     *h2o dens. *o3 dens.  ");
+	Output::End();
+
+	for(int i = 0; i < 34; i++)
 	{
-		Output::Begin();
-		Output::Repeat(10, ' ');
-		Output::Print(" atmospheric model identity :  ");
-		Output::End();
-
-		Output::Begin();
-		Output::Repeat(12, ' ');
-		ostringstream s1;
-		s1.setf(ios::fixed, ios::floatfield);
-		s1 << setprecision(3);
-		s1 << " user defined water content : uh2o=" << setw(9) << uw << " g/cm2 ";
-		Output::Print(s1.str());
-		Output::End();
-
-		Output::Begin();
-		Output::Repeat(12, ' ');
-		ostringstream s2;
-		s2.setf(ios::fixed, ios::floatfield);
-		s2 << setprecision(3);
-		s2 << " user defined ozone content : uo3 =" << setw(9) << uo3 << " cm-atm";
-		Output::Print(s2.str());
-		Output::End();
+	    Output::Begin();
+	    Output::Repeat(12, ' ');
+	    ostringstream s;
+	    s.setf(ios::fixed, ios::floatfield);
+	    s << setprecision(4);
+	    s << setw(9) << z[i] << "  ";
+	    s << setw(9) << p[i] << "  ";
+	    s << setw(9) << t[i] << "  ";
+	    s << setw(9) << wh[i] << "  ";
+	    s << setw(9) << wo[i] << "  ";
+	    s << ends;
+	    Output::Print(s.str());
+	    Output::End();
 	}
+    }
+    else 
+    {
+	Output::Begin();
+	Output::Repeat(10, ' ');
+	Output::Print(" atmospheric model identity :  ");
+	Output::End();
+
+	Output::Begin();
+	Output::Repeat(12, ' ');
+	ostringstream s1;
+	s1.setf(ios::fixed, ios::floatfield);
+	s1 << setprecision(3);
+	s1 << " user defined water content : uh2o=" << setw(9) << uw << " g/cm2 ";
+	Output::Print(s1.str());
+	Output::End();
+
+	Output::Begin();
+	Output::Repeat(12, ' ');
+	ostringstream s2;
+	s2.setf(ios::fixed, ios::floatfield);
+	s2 << setprecision(3);
+	s2 << " user defined ozone content : uo3 =" << setw(9) << uo3 << " cm-atm";
+	Output::Print(s2.str());
+	Output::End();
+    }
 
-	Output::Begin(); Output::End();
+    Output::Begin(); Output::End();
 }
 
 AtmosModel AtmosModel::Parse()
 {
-	AtmosModel atms;
-	atms.parse();
-	return atms;
+    AtmosModel atms;
+    atms.parse();
+    return atms;
 }

imagery/i.atcorr/GeomCond.cpp

@@ -36,23 +36,23 @@
 
 
 /*	To take into account the variation of the solar constant as a function 
-	of the Julian day. 
+  of the Julian day. 
 
-	return dsol		
-	dsol is a multiplicative factor to apply to the mean value of solar constant 
+  return dsol		
+  dsol is a multiplicative factor to apply to the mean value of solar constant 
 */
 float GeomCond::varsol ()
 {
 /* calculation of the variability of the solar constant during the year. 
    jday is the number of the day in the month   */
-  long int j;
-  if (month <= 2) j = (month - 1) * 31 + jday;
-  else if (month > 8) j = (month - 1) * 31 - (month - 2) / 2 - 2 + jday;
-  else j = (month - 1) * 31 - (month - 1) / 2 - 2 + jday;
+    long int j;
+    if (month <= 2) j = (month - 1) * 31 + jday;
+    else if (month > 8) j = (month - 1) * 31 - (month - 2) / 2 - 2 + jday;
+    else j = (month - 1) * 31 - (month - 1) / 2 - 2 + jday;
 
 /* Computing 2nd power */
-  double tmp = 1.f - cos ((float) (j - 4) * 0.9856f * M_PI / 180.f) * .01673f;
-  return 1.f / (float)(tmp * tmp);
+    double tmp = 1.f - cos ((float) (j - 4) * 0.9856f * M_PI / 180.f) * .01673f;
+    return 1.f / (float)(tmp * tmp);
 }
 
 
@@ -61,41 +61,41 @@ void GeomCond::landsat(float tu)
 {
 /*     warning !!! */
 /*     xlon and xlat are the coordinates of the scene center. */
-	avis = 0.f;
-	phiv = 0.f;
-	possol(tu);
+    avis = 0.f;
+    phiv = 0.f;
+    possol(tu);
 }
 
 /*
-To compute the solar azimuthal and zenithal angles (in degrees) for a point over
-the globe defined by its longitude and its latitude (in dec. degrees) for a day of the year (fixed by
-number of the month and number of the day in the month) at any Greenwich Meridian Time (GMT
-dec. hour).
+  To compute the solar azimuthal and zenithal angles (in degrees) for a point over
+  the globe defined by its longitude and its latitude (in dec. degrees) for a day of the year (fixed by
+  number of the month and number of the day in the month) at any Greenwich Meridian Time (GMT
+  dec. hour).
 */
 void GeomCond::possol(float tu)
 {
-	long int ia = 0;
-	long int nojour;
+    long int ia = 0;
+    long int nojour;
 /*     solar position (zenithal angle asol,azimuthal angle phi0 */
 /*                     in degrees) */
 /*     jday is the number of the day in the month */
-	day_number(ia, nojour);
-	pos_fft (nojour, tu);
-	if (asol > 90.f) fprintf(stderr, "The sun is not raised\n");
+    day_number(ia, nojour);
+    pos_fft (nojour, tu);
+    if (asol > 90.f) fprintf(stderr, "The sun is not raised\n");
 }
 
 void GeomCond::day_number(long int ia, long int& j)
 {
-	if (month <= 2)
-	{
-		j = (month - 1) * 31 + jday;
-		return;
+    if (month <= 2)
+    {
+	j = (month - 1) * 31 + jday;
+	return;
     }
 
-	if (month > 8) j = (month - 1) * 31 - (month - 2) / 2 - 2 + jday;
-	else j = (month - 1) * 31 - (month - 1) / 2 - 2 + jday;
+    if (month > 8) j = (month - 1) * 31 - (month - 2) / 2 - 2 + jday;
+    else j = (month - 1) * 31 - (month - 1) / 2 - 2 + jday;
 
-	if (ia != 0 && ia % 4 == 0) ++j;
+    if (ia != 0 && ia % 4 == 0) ++j;
 }
 
 /* returns the sign of the element */
@@ -103,112 +103,112 @@ void GeomCond::day_number(long int ia, long int& j)
 
 void GeomCond::pos_fft (long int j, float tu)
 {
-	/* Local variables */
-	double ah, et, az, caz, xla, tet, tsm, tsv, elev, azim, delta, amuzero;
+    /* Local variables */
+    double ah, et, az, caz, xla, tet, tsm, tsv, elev, azim, delta, amuzero;
 
-	/*     solar position (zenithal angle asol,azimuthal angle phi0 */
-	/*                     in degrees) */
-	/*     j is the day number in the year */
+    /*     solar position (zenithal angle asol,azimuthal angle phi0 */
+    /*                     in degrees) */
+    /*     j is the day number in the year */
 
-	/* mean solar time (heure decimale) */
-	tsm = tu + xlon / 15.;
-	xla = xlat * M_PI / 180.;
-	tet = (float)(j) * M_PI2 / 365.;
+    /* mean solar time (heure decimale) */
+    tsm = tu + xlon / 15.;
+    xla = xlat * M_PI / 180.;
+    tet = (float)(j) * M_PI2 / 365.;
 
-	/* time equation (in mn.dec) */
-	et = 7.5e-5f + 0.001868f * cos (tet) - 0.032077f * sin (tet) - 
-		 0.014615f * cos (tet * 2.f) - 0.040849f * sin (tet * 2.f);
+    /* time equation (in mn.dec) */
+    et = 7.5e-5f + 0.001868f * cos (tet) - 0.032077f * sin (tet) - 
+	0.014615f * cos (tet * 2.f) - 0.040849f * sin (tet * 2.f);
 
-	et = et * 12.f * 60.f / M_PI;
+    et = et * 12.f * 60.f / M_PI;
 
-	/* true solar time */
-	tsv = tsm + et / 60.f;
-	tsv += -12.f;
+    /* true solar time */
+    tsv = tsm + et / 60.f;
+    tsv += -12.f;
 
-	/* hour angle */
-	ah = tsv * 15.f * M_PI / 180.f;
+    /* hour angle */
+    ah = tsv * 15.f * M_PI / 180.f;
 
-	/* solar declination   (in radian) */
-	delta = 0.006918f - 0.399912f * cos (tet) + 0.070257f * sin (tet) - 
-			0.006758f * cos (tet * 2.f) + 9.07e-4f * sin (tet * 2.f) - 
-			0.002697f * cos (tet * 3.f) + 0.00148f * sin (tet * 3.f);
+    /* solar declination   (in radian) */
+    delta = 0.006918f - 0.399912f * cos (tet) + 0.070257f * sin (tet) - 
+	0.006758f * cos (tet * 2.f) + 9.07e-4f * sin (tet * 2.f) - 
+	0.002697f * cos (tet * 3.f) + 0.00148f * sin (tet * 3.f);
 
-	/* elevation,azimuth */
-	amuzero = sin (xla) * sin (delta) + cos (xla) * cos (delta) * cos (ah);
-	elev = asin (amuzero);
-	az = cos (delta) * sin (ah) / cos (elev);
+    /* elevation,azimuth */
+    amuzero = sin (xla) * sin (delta) + cos (xla) * cos (delta) * cos (ah);
+    elev = asin (amuzero);
+    az = cos (delta) * sin (ah) / cos (elev);
   
-	if (fabs (az) - 1.f > 0.f) az = SIGN(az);
+    if (fabs (az) - 1.f > 0.f) az = SIGN(az);
 
-	caz = (-cos (xla) * sin (delta) + sin (xla) * cos (delta) * cos (ah)) / cos (elev);
-	azim = asin (az);
-	if (caz <= 0.f) azim = M_PI - azim;
+    caz = (-cos (xla) * sin (delta) + sin (xla) * cos (delta) * cos (ah)) / cos (elev);
+    azim = asin (az);
+    if (caz <= 0.f) azim = M_PI - azim;
 
-	if (caz > 0.f && az <= 0.f) azim += M_PI2;
+    if (caz > 0.f && az <= 0.f) azim += M_PI2;
 
-	azim += M_PI;
-	if (azim > M_PI2) azim -= M_PI2;
+    azim += M_PI;
+    if (azim > M_PI2) azim -= M_PI2;
 	
-	elev = elev * 180. / M_PI;
+    elev = elev * 180. / M_PI;
 	
-	/*     conversion in degrees */
-	asol = (float)(90. - elev);
-	phi0 = (float)(azim * 180. / M_PI);
+    /*     conversion in degrees */
+    asol = (float)(90. - elev);
+    phi0 = (float)(azim * 180. / M_PI);
 }
 
 /*
-convert:
-1 = meteosat observation 
-2 = goes east observation
-3 = goes west observation
+  convert:
+  1 = meteosat observation 
+  2 = goes east observation
+  3 = goes west observation
 */
 void GeomCond::posobs(float tu, int nc, int nl)
 {
-	double yr, xr, alti;
-
-	if(igeom == 1) /* meteosat observation */
-	{
-		yr = nl - 1250.5;
-		xr = nc - 2500.5;
-		alti = 42164.0 - 6378.155;
-	} 
-	else if(igeom == 2) /* goes east observation */
-	{
-		yr = nl - 8665.5;
-		xr = nc - 6498.5;
-		alti = 42107.0 - 6378.155;
-	}
-	else /* goes west observation */
-	{
-      yr = nl - 8665.5;
-      xr = nc - 6498.5;
-      alti = 42147.0 - 6378.155;
-	}
-
-
-	const double re = 6378.155;
+    double yr, xr, alti;
+
+    if(igeom == 1) /* meteosat observation */
+    {
+	yr = nl - 1250.5;
+	xr = nc - 2500.5;
+	alti = 42164.0 - 6378.155;
+    } 
+    else if(igeom == 2) /* goes east observation */
+    {
+	yr = nl - 8665.5;
+	xr = nc - 6498.5;
+	alti = 42107.0 - 6378.155;
+    }
+    else /* goes west observation */
+    {
+	yr = nl - 8665.5;
+	xr = nc - 6498.5;
+	alti = 42147.0 - 6378.155;
+    }
+
+
+    const double re = 6378.155;
     const double aaa = 1. / 297.;
     const double rp = re / (1.f + aaa);
     const double cdr = M_PI / 180.;
     const double crd = 180. / M_PI;
 
-	double deltax;
-	double deltay;
-
-	if(igeom == 1) 
-	{
-		deltax = 18.0 / 5000.0;
-		deltay = 18.0 / 2500.0;
-	}
-	else
-	{
-		deltax = 18.0 / 12997.0;
-		deltay = 20.0 / 17331.0;
-	}
-
-	double x = xr * deltax * cdr;
+    double deltax;
+    double deltay;
+
+    if(igeom == 1) 
+    {
+	deltax = 18.0 / 5000.0;
+	deltay = 18.0 / 2500.0;
+    }
+    else
+    {
+	deltax = 18.0 / 12997.0;
+	deltay = 20.0 / 17331.0;
+    }
+
+    double x = xr * deltax * cdr;
     double y = yr * deltay * cdr;
-	double rs = re + alti;
+    double rs = re + alti;
     double tanx = tan(x);
     double tany = tan(y);
     double val1 = 1.0 + (tanx * tanx);
@@ -216,184 +216,184 @@ void GeomCond::posobs(float tu, int nc, int nl)
     double yk = rs / re;
     double cosx2 = 1. / (val1 * val2);
       
-	double sn, zt, xt, yt, teta, ylat, ylon;
-	if((1. / cosx2) > ((yk * yk) / (yk*yk - 1.)))
-	{
-		fprintf(stderr, "no possibility to compute lat. and long.\n");
-		return;
-	}
-	else
-	{
-      sn = (rs - (re * (sqrt((yk * yk) - (yk*yk - 1.) * (1. / cosx2))))) / (1. / cosx2);
-      zt = rs - sn;
-      xt = -(sn * tanx);
-      yt = sn * tany / cos(x);
-      teta = asin(yt / rp);
-      ylat = (atan(((tan(teta)) * rp) / re));
-      ylon = atan(xt / zt);
-	}
+    double sn, zt, xt, yt, teta, ylat, ylon;
+    if((1. / cosx2) > ((yk * yk) / (yk*yk - 1.)))
+    {
+	fprintf(stderr, "no possibility to compute lat. and long.\n");
+	return;
+    }
+    else
+    {
+	sn = (rs - (re * (sqrt((yk * yk) - (yk*yk - 1.) * (1. / cosx2))))) / (1. / cosx2);
+	zt = rs - sn;
+	xt = -(sn * tanx);
+	yt = sn * tany / cos(x);
+	teta = asin(yt / rp);
+	ylat = (atan(((tan(teta)) * rp) / re));
+	ylon = atan(xt / zt);
+    }
  
-	xlat = (float)(ylat * crd);
+    xlat = (float)(ylat * crd);
 
-	if(igeom == 1) xlon = (float)(ylon * crd);
-	else if(igeom == 2) xlon = (float)(ylon * crd - 75.);
-	else xlon = (float)(ylon * crd - 135.);
+    if(igeom == 1) xlon = (float)(ylon * crd);
+    else if(igeom == 2) xlon = (float)(ylon * crd - 75.);
+    else xlon = (float)(ylon * crd - 135.);
  
-	possol(tu);
+    possol(tu);
  
-	if(igeom == 1) ylon = xlon * M_PI / 180.;
-	else if(igeom == 2) ylon = xlon * M_PI / 180. + 75. * cdr;
-	else ylon = xlon * M_PI / 180. + 135. * cdr;
-
-	ylat = xlat * M_PI / 180.;
-	double gam = sqrt(((1. / cosx2) - 1.) * cosx2);
-	avis = (float)(asin((1. + alti / re) * (gam)) * 180. / M_PI);
-	phiv = (float)((atan2(tan(ylon),sin(ylat)) + M_PI) * 180. / M_PI);
+    if(igeom == 1) ylon = xlon * M_PI / 180.;
+    else if(igeom == 2) ylon = xlon * M_PI / 180. + 75. * cdr;
+    else ylon = xlon * M_PI / 180. + 135. * cdr;
+
+    ylat = xlat * M_PI / 180.;
+    double gam = sqrt(((1. / cosx2) - 1.) * cosx2);
+    avis = (float)(asin((1. + alti / re) * (gam)) * 180. / M_PI);
+    phiv = (float)((atan2(tan(ylon),sin(ylat)) + M_PI) * 180. / M_PI);
 }
 
 void GeomCond::posnoa(float tu, int nc, float xlonan, float campm, float hna)
 {
 /*     noaa 6 definition
-     orbite inclination ai in radians
-     hor mouv in rad/s  an
-     h/r=860/6378
-     campm allows the user to switch to pm platforms */
+       orbite inclination ai in radians
+       hor mouv in rad/s  an
+       h/r=860/6378
+       campm allows the user to switch to pm platforms */
  
-      const double r = 860. / 6378.155;
-      const double ai = 98.96 * M_PI / 180.;
-      const double an = 360. * M_PI / (6119. * 180.);
-      double ylonan = xlonan * M_PI / 180.;
-      double t = tu * 3600;
-      double hnam = hna;
-      hnam = hnam * 3600;
-      double u = t - hnam;
-      u = campm * u * an;
-      double delt = ((nc - (2048 + 1) / 2.) * 55.385 / ((2048. - 1) / 2.));
-      delt = campm * delt * M_PI / 180.;
-      avis = (float)asin((1 + r) * sin(delt));
-      double d = avis - delt;
-      double y = cos(d) * cos(ai) * sin(u) - sin(ai) * sin(d);
-      double z = cos(d) * sin(ai) * sin(u) + cos(ai) * sin(d);
-      double ylat = asin(z);
-      double cosy = cos(d) * cos(u) / cos(ylat);
-
-      double siny = y / cos(ylat);
-      double ylon = asin(siny);
-      if(cosy <= 0.)
-	  {
-		if(siny > 0) ylon = M_PI - ylon;
-		if(siny <= 0) ylon = -(M_PI + ylon);
-	  }
-      double ylo1 = ylon + ylonan - (t - hnam) * 2. * M_PI / 86400.;
-      xlat = (float)(ylat * 180. / M_PI);
-      xlon = (float)(ylo1 * 180. / M_PI);
+    const double r = 860. / 6378.155;
+    const double ai = 98.96 * M_PI / 180.;
+    const double an = 360. * M_PI / (6119. * 180.);
+    double ylonan = xlonan * M_PI / 180.;
+    double t = tu * 3600;
+    double hnam = hna;
+    hnam = hnam * 3600;
+    double u = t - hnam;
+    u = campm * u * an;
+    double delt = ((nc - (2048 + 1) / 2.) * 55.385 / ((2048. - 1) / 2.));
+    delt = campm * delt * M_PI / 180.;
+    avis = (float)asin((1 + r) * sin(delt));
+    double d = avis - delt;
+    double y = cos(d) * cos(ai) * sin(u) - sin(ai) * sin(d);
+    double z = cos(d) * sin(ai) * sin(u) + cos(ai) * sin(d);
+    double ylat = asin(z);
+    double cosy = cos(d) * cos(u) / cos(ylat);
+
+    double siny = y / cos(ylat);
+    double ylon = asin(siny);
+    if(cosy <= 0.)
+    {
+	if(siny > 0) ylon = M_PI - ylon;
+	if(siny <= 0) ylon = -(M_PI + ylon);
+    }
+    double ylo1 = ylon + ylonan - (t - hnam) * 2. * M_PI / 86400.;
+    xlat = (float)(ylat * 180. / M_PI);
+    xlon = (float)(ylo1 * 180. / M_PI);
  
 
 
-      possol(tu);
+    possol(tu);
  
-      double zlat = asin(sin(ai) * sin(u));
-      double zlon = atan2(cos(ai) * sin(u),cos(u));
-      if(nc != 1024)
-	  {
-		double xnum = sin(zlon - ylon) * cos(zlat) / sin(fabs(d));
-		double xden = (sin(zlat) - sin(ylat) * cos(d)) / cos(ylat) / sin(fabs(d));
-		phiv = (float)atan2(xnum,xden);
-	  }
-      else phiv = 0.;
-      phiv = (float)(phiv * 180. / M_PI);
-      avis = (float)(fabs(avis) * 180. / M_PI);
+    double zlat = asin(sin(ai) * sin(u));
+    double zlon = atan2(cos(ai) * sin(u),cos(u));
+    if(nc != 1024)
+    {
+	double xnum = sin(zlon - ylon) * cos(zlat) / sin(fabs(d));
+	double xden = (sin(zlat) - sin(ylat) * cos(d)) / cos(ylat) / sin(fabs(d));
+	phiv = (float)atan2(xnum,xden);
+    }
+    else phiv = 0.;
+    phiv = (float)(phiv * 180. / M_PI);
+    avis = (float)(fabs(avis) * 180. / M_PI);
 }
 
 void GeomCond::parse()
 {
-	cin >> igeom;
+    cin >> igeom;
+    cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
+
+    float campm = -1.0f;	/* initialize in case igeom == 5 */
+    float tu, xlonan, hna;
+    int nc, nl;
+
+    switch(igeom)
+    {
+    case 0: /* internal format */
+    {
+	cin >> asol;
+	cin >> phi0;
+	cin >> avis;
+	cin >> phiv;
+	cin >> month;
+	cin >> jday;
 	cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
+	break;
+    }
+    case 1: /* meteosat observation */
+    case 2: /* goes east observation */
+    case 3: /* goes west observation  */
+    {
+	cin >> month;
+	cin >> jday;
+	cin >> tu;
+	cin >> nc;
+	cin >> nl;
+	cin.ignore(numeric_limits<int>::max(),'\n');
+	posobs(tu, nc, nl);
+	break;
+    }
+    case 4: campm = 1.0f;
+    case 5: 
+    {
+	cin >> month;
+	cin >> jday;
+	cin >> tu;
+	cin >> nc;
+	cin >> xlonan;
+	cin >> hna;
+	cin.ignore(numeric_limits<int>::max(),'\n');
+	posnoa(tu, nc, xlonan, campm, hna);
+	break;
+    }
+    case 6: /* hrv   ( spot )    * enter month,day,hh.ddd,long.,lat. */
+    case 7: /* tm    ( landsat ) * enter month,day,hh.ddd,long.,lat. */
+    case 8: /* etm+  ( landsat7) * enter month,day,hh.ddd,long.,lat. */
+    {
+	cin >> month;
+	cin >> jday;
+	cin >> tu;
+	cin >> xlon;
+	cin >> xlat;
+	cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
+	landsat(tu);
+	break;
+    }
+    default: fprintf(stderr, "Unsupported format.\n");
+    }
+
 
-	float campm = -1.0f;	/* initialize in case igeom == 5 */
-	float tu, xlonan, hna;
-	int nc, nl;
-
-	switch(igeom)
-	{
-	case 0: /* internal format */
-		{
-			cin >> asol;
-			cin >> phi0;
-			cin >> avis;
-			cin >> phiv;
-			cin >> month;
-			cin >> jday;
-			cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
-			break;
-		}
-	case 1: /* meteosat observation */
-	case 2: /* goes east observation */
-	case 3: /* goes west observation  */
-		{
-			cin >> month;
-			cin >> jday;
-			cin >> tu;
-			cin >> nc;
-			cin >> nl;
-			cin.ignore(numeric_limits<int>::max(),'\n');
-			posobs(tu, nc, nl);
-			break;
-		}
-	case 4: campm = 1.0f;
-	case 5: 
-		{
-			cin >> month;
-			cin >> jday;
-			cin >> tu;
-			cin >> nc;
-			cin >> xlonan;
-			cin >> hna;
-			cin.ignore(numeric_limits<int>::max(),'\n');
-			posnoa(tu, nc, xlonan, campm, hna);
-			break;
-		}
-	case 6: /* hrv   ( spot )    * enter month,day,hh.ddd,long.,lat. */
-	case 7: /* tm    ( landsat ) * enter month,day,hh.ddd,long.,lat. */
-	case 8: /* etm+  ( landsat7) * enter month,day,hh.ddd,long.,lat. */
-		{
-			cin >> month;
-			cin >> jday;
-			cin >> tu;
-			cin >> xlon;
-			cin >> xlat;
-			cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
-			landsat(tu);
-			break;
-		}
-	default: fprintf(stderr, "Unsupported format.\n");
-	}
-
-
-	/* ********************************************************************** */
-	/*                                                                        */
-	/*                                 / scattered direction                  */
-	/*                               /                                        */
-	/*                             /                                          */
-	/*                           / adif                                       */
-	/*    incident   + + + + + + + + + + + + + + +                            */
-	/*    direction                                                           */
-	/*                                                                        */
-	/* ********************************************************************** */
-	phi = (float)fabs(phiv - phi0);
-	phirad = (phi0 - phiv) * (float)M_PI / 180.f;
-	if (phirad < 0.f) phirad += (float)M_PI2;
-	if (phirad > M_PI2) phirad -= (float)M_PI2;
-
-	xmus = (float)cos (asol * M_PI / 180.f);
-	xmuv = (float)cos (avis * M_PI / 180.f);
-	xmup = (float)cos (phirad);
-	xmud = -xmus * xmuv - (float)sqrt (1.f - xmus * xmus) * (float)sqrt (1.f - xmuv * xmuv) * xmup;
-
-	/* test vermote bug */
-	if (xmud > 1.f)  xmud = 1.f;
-	if (xmud < -1.f) xmud = -1.f;
-	adif = (float)acos (xmud) * 180.f / (float)M_PI;
+    /* ********************************************************************** */
+    /*                                                                        */
+    /*                                 / scattered direction                  */
+    /*                               /                                        */
+    /*                             /                                          */
+    /*                           / adif                                       */
+    /*    incident   + + + + + + + + + + + + + + +                            */
+    /*    direction                                                           */
+    /*                                                                        */
+    /* ********************************************************************** */
+    phi = (float)fabs(phiv - phi0);
+    phirad = (phi0 - phiv) * (float)M_PI / 180.f;
+    if (phirad < 0.f) phirad += (float)M_PI2;
+    if (phirad > M_PI2) phirad -= (float)M_PI2;
+
+    xmus = (float)cos (asol * M_PI / 180.f);
+    xmuv = (float)cos (avis * M_PI / 180.f);
+    xmup = (float)cos (phirad);
+    xmud = -xmus * xmuv - (float)sqrt (1.f - xmus * xmus) * (float)sqrt (1.f - xmuv * xmuv) * xmup;
+
+    /* test vermote bug */
+    if (xmud > 1.f)  xmud = 1.f;
+    if (xmud < -1.f) xmud = -1.f;
+    adif = (float)acos (xmud) * 180.f / (float)M_PI;
 
     dsol = varsol();
 }
@@ -401,75 +401,75 @@ void GeomCond::parse()
 /* ---- print geometrical conditions ---- */
 void GeomCond::print()
 {
-	static const string etiq1[9] = {
-     string(" user defined conditions     "),
-     string(" meteosat observation        "),
-     string(" goes east observation       "),
-     string(" goes west observation       "),
-     string(" avhrr (AM noaa) observation "),
-     string(" avhrr (PM noaa) observation "),
-     string(" h.r.v.   observation        "),
-     string(" t.m.     observation        "),
-     string(" etm+     observation        ")
-	};
-
-	static const string head(" geometrical conditions identity  ");
-	static const string line(" -------------------------------  ");
-	Output::Begin(); Output::Repeat(22,' '); Output::Print(head); Output::End();
-	Output::Begin(); Output::Repeat(22,' '); Output::Print(line); Output::End();
+    static const string etiq1[9] = {
+	string(" user defined conditions     "),
+	string(" meteosat observation        "),
+	string(" goes east observation       "),
+	string(" goes west observation       "),
+	string(" avhrr (AM noaa) observation "),
+	string(" avhrr (PM noaa) observation "),
+	string(" h.r.v.   observation        "),
+	string(" t.m.     observation        "),
+	string(" etm+     observation        ")
+    };
+
+    static const string head(" geometrical conditions identity  ");
+    static const string line(" -------------------------------  ");
+    Output::Begin(); Output::Repeat(22,' '); Output::Print(head); Output::End();
+    Output::Begin(); Output::Repeat(22,' '); Output::Print(line); Output::End();
 
 	
-	Output::Begin(); Output::Repeat(22,' '); Output::Print(etiq1[igeom]); Output::End();
-	Output::Begin(); Output::End();
+    Output::Begin(); Output::Repeat(22,' '); Output::Print(etiq1[igeom]); Output::End();
+    Output::Begin(); Output::End();
 
 	
-	Output::Begin(); Output::Repeat(2,' ');
-	ostringstream s1;
-	s1.setf(ios::fixed, ios::floatfield);
-	s1 << " month: " << month << " day: " << jday;
-	s1 << ends;
-	Output::Print(s1.str());
-	Output::End();
+    Output::Begin(); Output::Repeat(2,' ');
+    ostringstream s1;
+    s1.setf(ios::fixed, ios::floatfield);
+    s1 << " month: " << month << " day: " << jday;
+    s1 << ends;
+    Output::Print(s1.str());
+    Output::End();
 
 
-	Output::Begin(); Output::Repeat(2,' ');
-	ostringstream s2;
-	s2.setf(ios::fixed, ios::floatfield);
-	s2 << setprecision(2);
+    Output::Begin(); Output::Repeat(2,' ');
+    ostringstream s2;
+    s2.setf(ios::fixed, ios::floatfield);
+    s2 << setprecision(2);
 
 
-	s2 << " solar zenith angle:  " << setw(6) << asol << " deg ";
-	s2 << " solar azimuthal angle:      " << setw(6) << phi0 << " deg";
-	s2 << ends;
-	Output::Print(s2.str());
-	Output::End();
+    s2 << " solar zenith angle:  " << setw(6) << asol << " deg ";
+    s2 << " solar azimuthal angle:      " << setw(6) << phi0 << " deg";
+    s2 << ends;
+    Output::Print(s2.str());
+    Output::End();
 
 	
-	Output::Begin(); Output::Repeat(2,' ');
-	ostringstream s3;
-	s3.setf(ios::fixed, ios::floatfield);
-	s3 << setprecision(2);
-	s3 << " view zenith angle:   " << setw(6) << avis << " deg ";
-	s3 << " view azimuthal angle:       " << setw(6) << phiv << " deg ";
-	s3 << ends;
-	Output::Print(s3.str());
-	Output::End();
-	Output::Begin(); Output::Repeat(2,' ');
-	ostringstream s4;
-	s4.setf(ios::fixed, ios::floatfield);
-	s4 << setprecision(2);
-	s4 << " scattering angle:    " << setw(6) << adif << " deg ";
-	s4 << " azimuthal angle difference: " << setw(6) << phi << " deg ";
-	s4 << ends;
-	Output::Print(s4.str());
-	Output::End();
+    Output::Begin(); Output::Repeat(2,' ');
+    ostringstream s3;
+    s3.setf(ios::fixed, ios::floatfield);
+    s3 << setprecision(2);
+    s3 << " view zenith angle:   " << setw(6) << avis << " deg ";
+    s3 << " view azimuthal angle:       " << setw(6) << phiv << " deg ";
+    s3 << ends;
+    Output::Print(s3.str());
+    Output::End();
+    Output::Begin(); Output::Repeat(2,' ');
+    ostringstream s4;
+    s4.setf(ios::fixed, ios::floatfield);
+    s4 << setprecision(2);
+    s4 << " scattering angle:    " << setw(6) << adif << " deg ";
+    s4 << " azimuthal angle difference: " << setw(6) << phi << " deg ";
+    s4 << ends;
+    Output::Print(s4.str());
+    Output::End();
 	
-	Output::Begin(); Output::End();
+    Output::Begin(); Output::End();
 }
 
 GeomCond GeomCond::Parse()
 {
-	GeomCond geom;
-	geom.parse();
-	return geom;
+    GeomCond geom;
+    geom.parse();
+    return geom;
 }

imagery/i.atcorr/Interp.cpp

@@ -2,202 +2,202 @@
 #include "Interp.h"
 
 void interp (const int iaer, const int idatmp, 
-			 const float wl, const float taer55, 
-			 const float taer55p, const float xmud, 
-			 InterpStruct& is)
+	     const float wl, const float taer55, 
+	     const float taer55p, const float xmud, 
+	     InterpStruct& is)
 {
 /*     that for the atmosphere :
-     the reflectances
-                     rayleigh                             = rorayl
-                     aerosols                             = roaero
-                     mixing                               = romix
-     the downward transmittances
-                     rayleigh                             = dtotr
-                     aerosols                             = dtota
-                     total                                = dtott
-     the upward transmittances
-                     rayleigh                             = utotr
-                     aerosols                             = utota
-                     total                                = utott
-     the spherical albedos
-                     rayleigh                             = asray
-                     aerosols                             = asaer
-                     total                                = astot
-     the optical thickness of total atmosphere
-                     rayleigh                             = tray
-                     aerosols                             = taer
-     the optical thickness of the atmosphere above the plane
-                     rayleigh                             = is.trayp
-                     aerosols                             = taerp
-     the tsca of the aerosols (god dammed it)
-                     total atmosphere                     = tsca */
+       the reflectances
+       rayleigh                             = rorayl
+       aerosols                             = roaero
+       mixing                               = romix
+       the downward transmittances
+       rayleigh                             = dtotr
+       aerosols                             = dtota
+       total                                = dtott
+       the upward transmittances
+       rayleigh                             = utotr
+       aerosols                             = utota
+       total                                = utott
+       the spherical albedos
+       rayleigh                             = asray
+       aerosols                             = asaer
+       total                                = astot
+       the optical thickness of total atmosphere
+       rayleigh                             = tray
+       aerosols                             = taer
+       the optical thickness of the atmosphere above the plane
+       rayleigh                             = is.trayp
+       aerosols                             = taerp
+       the tsca of the aerosols (god dammed it)
+       total atmosphere                     = tsca */
       
-	int linf = 0;
-	for(int i = 0; i < 9; i++) if(wl > sixs_disc.wldis[i] && wl <= sixs_disc.wldis[i+1]) linf = i;
+    int linf = 0;
+    for(int i = 0; i < 9; i++) if(wl > sixs_disc.wldis[i] && wl <= sixs_disc.wldis[i+1]) linf = i;
 	
     if(wl > sixs_disc.wldis[9]) linf = 8;
-	int lsup = linf + 1;
+    int lsup = linf + 1;
 
 
-	/*    interpolation in function of wavelength for scattering
-	    atmospheric functions from discrete values at sixs_disc.wldis */
+    /*    interpolation in function of wavelength for scattering
+	  atmospheric functions from discrete values at sixs_disc.wldis */
  
-	float alphaa = 0;
-	float betaa = 0;
-	float alphar = 0;
-	float betar = 0;
-	float alphac = 0;
-	float betac = 0;
-	is.phaa = 0;
-	is.roaero = 0;
-	is.dtota = 1;
-	is.utota = 1;
-	is.asaer = 0;
-	is.taer = 0;
-	is.taerp = 0;
-	float coef = (float)log(sixs_disc.wldis[lsup] / sixs_disc.wldis[linf]);
+    float alphaa = 0;
+    float betaa = 0;
+    float alphar = 0;
+    float betar = 0;
+    float alphac = 0;
+    float betac = 0;
+    is.phaa = 0;
+    is.roaero = 0;
+    is.dtota = 1;
+    is.utota = 1;
+    is.asaer = 0;
+    is.taer = 0;
+    is.taerp = 0;
+    float coef = (float)log(sixs_disc.wldis[lsup] / sixs_disc.wldis[linf]);
     float wlinf = sixs_disc.wldis[linf];
 
     if(iaer != 0)
-	{
-		alphaa = (float)(log(sixs_aer.phase[lsup] / sixs_aer.phase[linf]) / coef);
-		betaa = (float)(sixs_aer.phase[linf] / pow(wlinf,alphaa));
-		is.phaa = (float)(betaa * pow(wl,alphaa));
-	}
+    {
+	alphaa = (float)(log(sixs_aer.phase[lsup] / sixs_aer.phase[linf]) / coef);
+	betaa = (float)(sixs_aer.phase[linf] / pow(wlinf,alphaa));
+	is.phaa = (float)(betaa * pow(wl,alphaa));
+    }
 
-	float d2 = 2 + delta;
+    float d2 = 2 + delta;
     is.phar = (2 * (1 - delta) / d2) * .75f * (1 + xmud * xmud) + 3 * delta / d2;
-	if(idatmp == 0)
-	{
-		betar = 0;
+    if(idatmp == 0)
+    {
+	betar = 0;
         betaa = 0;
         betac = 0;
+    }
+    else
+    {
+	if(sixs_disc.roatm[0][linf] < 0.001)
+	{
+	    is.rorayl = sixs_disc.roatm[0][linf] + (sixs_disc.roatm[0][lsup] - sixs_disc.roatm[0][linf])
+		* (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
 	}
 	else
 	{
-		if(sixs_disc.roatm[0][linf] < 0.001)
-		{
-			is.rorayl = sixs_disc.roatm[0][linf] + (sixs_disc.roatm[0][lsup] - sixs_disc.roatm[0][linf])
-				   * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
-		}
-		else
-		{
-			alphar = (float)(log(sixs_disc.roatm[0][lsup] / sixs_disc.roatm[0][linf]) / coef);
-			betar = (float)(sixs_disc.roatm[0][linf] / pow(wlinf,alphar));
-			is.rorayl = (float)(betar * pow(wl,alphar));
-		}
-
-		if(sixs_disc.roatm[1][linf] < 0.001)
-		{
-			is.romix = sixs_disc.roatm[1][linf] + (sixs_disc.roatm[1][lsup] - sixs_disc.roatm[1][linf])
-				  * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
-		}
-        else
-		{
-			alphac = (float)(log(sixs_disc.roatm[1][lsup] / sixs_disc.roatm[1][linf]) / coef);
-			betac = (float)(sixs_disc.roatm[1][linf] / pow(wlinf,alphac));
-			is.romix = (float)(betac * pow(wl,alphac));
-		}
-
-		if(iaer != 0)
-		{
-
-			if(sixs_disc.roatm[2][linf] < 0.001)
-			{
-				is.roaero = sixs_disc.roatm[2][linf]+(sixs_disc.roatm[2][lsup] - sixs_disc.roatm[2][linf])
-				       * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
-			}
-			else
-			{
-				alphaa = (float)(log(sixs_disc.roatm[2][lsup] / sixs_disc.roatm[2][linf]) / coef);
-				betaa = (float)(sixs_disc.roatm[2][linf] / pow(wlinf,alphaa));
-				is.roaero = (float)(betaa * pow(wl,alphaa));
-			}
-		}
+	    alphar = (float)(log(sixs_disc.roatm[0][lsup] / sixs_disc.roatm[0][linf]) / coef);
+	    betar = (float)(sixs_disc.roatm[0][linf] / pow(wlinf,alphar));
+	    is.rorayl = (float)(betar * pow(wl,alphar));
 	}
 
-	alphar = (float)(log(sixs_disc.trayl[lsup] / sixs_disc.trayl[linf]) / coef);
-    betar = (float)(sixs_disc.trayl[linf] / pow(wlinf,alphar));
-    is.tray = (float)(betar * pow(wl,alphar));
-    
-	if (idatmp != 0)
+	if(sixs_disc.roatm[1][linf] < 0.001)
 	{
-		alphar = (float)(log(sixs_disc.traypl[lsup] / sixs_disc.traypl[linf]) / coef);
-        betar = (float)(sixs_disc.traypl[linf] / pow(wlinf,alphar));
-        is.trayp = (float)(betar * pow(wl,alphar));
+	    is.romix = sixs_disc.roatm[1][linf] + (sixs_disc.roatm[1][lsup] - sixs_disc.roatm[1][linf])
+		* (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
 	}
-    else is.trayp = 0;
-
-	if(iaer != 0)
+        else
 	{
-		alphaa = (float)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] / (sixs_aer.ext[linf] * sixs_aer.ome[linf])) / coef);
-		betaa = (float)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
-		is.tsca = (float)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
-		alphaa = (float)(log(sixs_aer.ext[lsup] / sixs_aer.ext[linf]) / coef);
-		betaa = (float)(sixs_aer.ext[linf] / pow(wlinf,alphaa));
-		is.taerp = (float)(taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
-		is.taer = (float)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+	    alphac = (float)(log(sixs_disc.roatm[1][lsup] / sixs_disc.roatm[1][linf]) / coef);
+	    betac = (float)(sixs_disc.roatm[1][linf] / pow(wlinf,alphac));
+	    is.romix = (float)(betac * pow(wl,alphac));
 	}
 
-	float drinf = sixs_disc.dtdif[0][linf] + sixs_disc.dtdir[0][linf];
-	float drsup = sixs_disc.dtdif[0][lsup] + sixs_disc.dtdir[0][lsup];
-	alphar = (float)(log(drsup / drinf) / coef);
-	betar = (float)(drinf / pow(wlinf,alphar));
-	is.dtotr = (float)(betar * pow(wl,alphar));
-	float dtinf = sixs_disc.dtdif[1][linf] + sixs_disc.dtdir[1][linf];
-	float dtsup = sixs_disc.dtdif[1][lsup] + sixs_disc.dtdir[1][lsup];
-	alphac = (float)(log((dtsup * drinf) / (dtinf * drsup)) / coef);
-	betac = (float)((dtinf / drinf) / pow(wlinf,alphac));
-	float dtotc = (float)(betac * pow(wl,alphac));
-	float dainf = sixs_disc.dtdif[2][linf] + sixs_disc.dtdir[2][linf];
-	float dasup = sixs_disc.dtdif[2][lsup] + sixs_disc.dtdir[2][lsup];
-
-	if(iaer != 0) 
+	if(iaer != 0)
 	{
-		alphaa = (float)(log(dasup / dainf) / coef);
-		betaa = (float)(dainf / pow(wlinf,alphaa));
-		is.dtota = (float)(betaa * pow(wl,alphaa));
+
+	    if(sixs_disc.roatm[2][linf] < 0.001)
+	    {
+		is.roaero = sixs_disc.roatm[2][linf]+(sixs_disc.roatm[2][lsup] - sixs_disc.roatm[2][linf])
+		    * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
+	    }
+	    else
+	    {
+		alphaa = (float)(log(sixs_disc.roatm[2][lsup] / sixs_disc.roatm[2][linf]) / coef);
+		betaa = (float)(sixs_disc.roatm[2][linf] / pow(wlinf,alphaa));
+		is.roaero = (float)(betaa * pow(wl,alphaa));
+	    }
 	}
+    }
 
-	is.dtott = dtotc * is.dtotr;
-	float urinf = sixs_disc.utdif[0][linf] + sixs_disc.utdir[0][linf];
-	float ursup = sixs_disc.utdif[0][lsup] + sixs_disc.utdir[0][lsup];
-	alphar = (float)(log(ursup / urinf) / coef);
-	betar = (float)(urinf / pow(wlinf,alphar));
-	is.utotr = (float)(betar * pow(wl,alphar));
-	float utinf = sixs_disc.utdif[1][linf] + sixs_disc.utdir[1][linf];
-	float utsup = sixs_disc.utdif[1][lsup] + sixs_disc.utdir[1][lsup];
-	alphac = (float)(log((utsup * urinf) / (utinf * ursup)) / coef);
-	betac = (float)((utinf / urinf) / pow(wlinf,alphac));
-	float utotc = (float)(betac * pow(wl,alphac));
-	float uainf = sixs_disc.utdif[2][linf] + sixs_disc.utdir[2][linf];
-	float uasup = sixs_disc.utdif[2][lsup] + sixs_disc.utdir[2][lsup];
-	is.utott = utotc * is.utotr;
+    alphar = (float)(log(sixs_disc.trayl[lsup] / sixs_disc.trayl[linf]) / coef);
+    betar = (float)(sixs_disc.trayl[linf] / pow(wlinf,alphar));
+    is.tray = (float)(betar * pow(wl,alphar));
+    
+    if (idatmp != 0)
+    {
+	alphar = (float)(log(sixs_disc.traypl[lsup] / sixs_disc.traypl[linf]) / coef);
+        betar = (float)(sixs_disc.traypl[linf] / pow(wlinf,alphar));
+        is.trayp = (float)(betar * pow(wl,alphar));
+    }
+    else is.trayp = 0;
 
-	if(iaer != 0)
-	{
-		alphaa = (float)(log(uasup / uainf) / coef);
-		betaa = (float)(uainf / pow(wlinf,alphaa));
-		is.utota = (float)(betaa * pow(wl,alphaa));
-	}
+    if(iaer != 0)
+    {
+	alphaa = (float)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] / (sixs_aer.ext[linf] * sixs_aer.ome[linf])) / coef);
+	betaa = (float)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
+	is.tsca = (float)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+	alphaa = (float)(log(sixs_aer.ext[lsup] / sixs_aer.ext[linf]) / coef);
+	betaa = (float)(sixs_aer.ext[linf] / pow(wlinf,alphaa));
+	is.taerp = (float)(taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+	is.taer = (float)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+    }
+
+    float drinf = sixs_disc.dtdif[0][linf] + sixs_disc.dtdir[0][linf];
+    float drsup = sixs_disc.dtdif[0][lsup] + sixs_disc.dtdir[0][lsup];
+    alphar = (float)(log(drsup / drinf) / coef);
+    betar = (float)(drinf / pow(wlinf,alphar));
+    is.dtotr = (float)(betar * pow(wl,alphar));
+    float dtinf = sixs_disc.dtdif[1][linf] + sixs_disc.dtdir[1][linf];
+    float dtsup = sixs_disc.dtdif[1][lsup] + sixs_disc.dtdir[1][lsup];
+    alphac = (float)(log((dtsup * drinf) / (dtinf * drsup)) / coef);
+    betac = (float)((dtinf / drinf) / pow(wlinf,alphac));
+    float dtotc = (float)(betac * pow(wl,alphac));
+    float dainf = sixs_disc.dtdif[2][linf] + sixs_disc.dtdir[2][linf];
+    float dasup = sixs_disc.dtdif[2][lsup] + sixs_disc.dtdir[2][lsup];
+
+    if(iaer != 0) 
+    {
+	alphaa = (float)(log(dasup / dainf) / coef);
+	betaa = (float)(dainf / pow(wlinf,alphaa));
+	is.dtota = (float)(betaa * pow(wl,alphaa));
+    }
+
+    is.dtott = dtotc * is.dtotr;
+    float urinf = sixs_disc.utdif[0][linf] + sixs_disc.utdir[0][linf];
+    float ursup = sixs_disc.utdif[0][lsup] + sixs_disc.utdir[0][lsup];
+    alphar = (float)(log(ursup / urinf) / coef);
+    betar = (float)(urinf / pow(wlinf,alphar));
+    is.utotr = (float)(betar * pow(wl,alphar));
+    float utinf = sixs_disc.utdif[1][linf] + sixs_disc.utdir[1][linf];
+    float utsup = sixs_disc.utdif[1][lsup] + sixs_disc.utdir[1][lsup];
+    alphac = (float)(log((utsup * urinf) / (utinf * ursup)) / coef);
+    betac = (float)((utinf / urinf) / pow(wlinf,alphac));
+    float utotc = (float)(betac * pow(wl,alphac));
+    float uainf = sixs_disc.utdif[2][linf] + sixs_disc.utdir[2][linf];
+    float uasup = sixs_disc.utdif[2][lsup] + sixs_disc.utdir[2][lsup];
+    is.utott = utotc * is.utotr;
 
-	float arinf = sixs_disc.sphal[0][linf];
-	float arsup = sixs_disc.sphal[0][lsup];
-	alphar = (float)(log(arsup / arinf) / coef);
-	betar = (float)(arinf / pow(wlinf,alphar));
-	is.asray = (float)(betar * pow(wl,alphar));
-	float atinf = sixs_disc.sphal[1][linf];
-	float atsup = sixs_disc.sphal[1][lsup];
-	alphac = (float)(log(atsup / atinf) / coef);
-	betac = (float)(atinf / pow(wlinf,alphac));
-	is.astot = (float)(betac * pow(wl,alphac));
-	float aainf = sixs_disc.sphal[2][linf];
-	float aasup = sixs_disc.sphal[2][lsup];
+    if(iaer != 0)
+    {
+	alphaa = (float)(log(uasup / uainf) / coef);
+	betaa = (float)(uainf / pow(wlinf,alphaa));
+	is.utota = (float)(betaa * pow(wl,alphaa));
+    }
+
+    float arinf = sixs_disc.sphal[0][linf];
+    float arsup = sixs_disc.sphal[0][lsup];
+    alphar = (float)(log(arsup / arinf) / coef);
+    betar = (float)(arinf / pow(wlinf,alphar));
+    is.asray = (float)(betar * pow(wl,alphar));
+    float atinf = sixs_disc.sphal[1][linf];
+    float atsup = sixs_disc.sphal[1][lsup];
+    alphac = (float)(log(atsup / atinf) / coef);
+    betac = (float)(atinf / pow(wlinf,alphac));
+    is.astot = (float)(betac * pow(wl,alphac));
+    float aainf = sixs_disc.sphal[2][linf];
+    float aasup = sixs_disc.sphal[2][lsup];
 
-	if(iaer != 0)
-	{
-		alphaa = (float)(log(aasup / aainf) / coef);
-		betaa = (float)(aainf / pow(wlinf,alphaa));
-		is.asaer = (float)(betaa * pow(wl,alphaa));
-	}
+    if(iaer != 0)
+    {
+	alphaa = (float)(log(aasup / aainf) / coef);
+	betaa = (float)(aainf / pow(wlinf,alphaa));
+	is.asaer = (float)(betaa * pow(wl,alphaa));
+    }
 }

imagery/i.atcorr/Transform.cpp

@@ -3,135 +3,135 @@
 
 void EtmDN(int iwave, float asol, bool before, float &lmin, float &lmax)
 {
-	if (before)		/* ETM+ digital numbers taken before July 1, 2000 */
-	{
-		switch(iwave)
-		{
-		case 61:
-			{
-				lmin = -6.2f;
-				lmax = 194.3f;
-				break;
-			}
-
-		case 62:
-			{
-				lmin = -6.0f;
-				lmax = 202.4f;
-				break;
-			}
-
-		case 63:
-			{
-				lmin = -4.5f;
-				lmax = 158.6f;
-				break;
-			}
-
-		case 64:
-			{
-				if (asol < 45.)
-				{
-					lmin = -4.5f;
-					lmax = 235.0f;
-				}
-				else
-				{
-					lmin = -4.5f;
-					lmax = 157.5f;
-				}
-				break;
-			}
-
-		case 65:
-			{
-				lmin = -1.0f;
-				lmax = 31.76f;
-				break;
-			}
-
-		case 66:
-			{
-				lmin = -0.35f;
-				lmax = 10.932f;
-				break;
-			}
-
-		case 67:
-			{
-				lmin = -5.0f;
-				lmax = 244.00f;
-				break;
-			}
-		}
-	}
-	else		/* ETM+ digital numbers taken after July 1, 2000 */
-	{
-		switch(iwave)
-		{
-		case 61:
-			{
-				lmin = -6.2f;
-				lmax = 191.6f;
-				break;
-			}
-
-		case 62:
-			{
-				lmin = -6.4f;
-				lmax = 196.5f;
-				break;
-			}
-
-		case 63:
-			{
-				lmin = -5.0f;
-				lmax = 152.9f;
-				break;
-			}
-
-		case 64:
-			{
-				if (asol < 45.)
-				{
-					lmin = -5.1f;
-					lmax = 241.1f;
-				}
-				else
-				{
-					lmin = -5.1f;
-					lmax = 157.4f;
-				}
-				break;
-			}
-
-		case 65:
-			{
-				lmin = -1.0f;
-				lmax = 31.06f;
-				break;
-			}
-
-		case 66:
-			{
-				lmin = -0.35f;
-				lmax = 10.80f;
-				break;
-			}
-
-		case 67:
-			{
-				lmin = -4.7f;
-				lmax = 243.1f;
-				break;
-			}
-		}
+    if (before)		/* ETM+ digital numbers taken before July 1, 2000 */
+    {
+	switch(iwave)
+	{
+	case 61:
+	{
+	    lmin = -6.2f;
+	    lmax = 194.3f;
+	    break;
+	}
+
+	case 62:
+	{
+	    lmin = -6.0f;
+	    lmax = 202.4f;
+	    break;
+	}
+
+	case 63:
+	{
+	    lmin = -4.5f;
+	    lmax = 158.6f;
+	    break;
+	}
+
+	case 64:
+	{
+	    if (asol < 45.)
+	    {
+		lmin = -4.5f;
+		lmax = 235.0f;
+	    }
+	    else
+	    {
+		lmin = -4.5f;
+		lmax = 157.5f;
+	    }
+	    break;
+	}
+
+	case 65:
+	{
+	    lmin = -1.0f;
+	    lmax = 31.76f;
+	    break;
+	}
+
+	case 66:
+	{
+	    lmin = -0.35f;
+	    lmax = 10.932f;
+	    break;
+	}
+
+	case 67:
+	{
+	    lmin = -5.0f;
+	    lmax = 244.00f;
+	    break;
+	}
+	}
+    }
+    else		/* ETM+ digital numbers taken after July 1, 2000 */
+    {
+	switch(iwave)
+	{
+	case 61:
+	{
+	    lmin = -6.2f;
+	    lmax = 191.6f;
+	    break;
 	}
+
+	case 62:
+	{
+	    lmin = -6.4f;
+	    lmax = 196.5f;
+	    break;
+	}
+
+	case 63:
+	{
+	    lmin = -5.0f;
+	    lmax = 152.9f;
+	    break;
+	}
+
+	case 64:
+	{
+	    if (asol < 45.)
+	    {
+		lmin = -5.1f;
+		lmax = 241.1f;
+	    }
+	    else
+	    {
+		lmin = -5.1f;
+		lmax = 157.4f;
+	    }
+	    break;
+	}
+
+	case 65:
+	{
+	    lmin = -1.0f;
+	    lmax = 31.06f;
+	    break;
+	}
+
+	case 66:
+	{
+	    lmin = -0.35f;
+	    lmax = 10.80f;
+	    break;
+	}
+
+	case 67:
+	{
+	    lmin = -4.7f;
+	    lmax = 243.1f;
+	    break;
+	}
+	}
+    }
 }
 
 /* Assuming input value between 0 and 1
-if rad is true, idn should first be converted to a reflectance value
-returns adjusted value also between 0 and 1 */
+   if rad is true, idn should first be converted to a reflectance value
+   returns adjusted value also between 0 and 1 */
 float transform(const TransformInput ti, InputMask imask, float idn)
 {
     /* convert from radiance to reflectance */
@@ -148,40 +148,40 @@ float transform(const TransformInput ti, InputMask imask, float idn)
     }
     if(imask & RADIANCE) idn += (float)M_PI * idn * 255.f * ti.sb / ti.xmus / ti.seb;
           
-	float rapp = idn;
+    float rapp = idn;
     float ainrpix = ti.ainr[0][0];
-	float xa = 0.0f;
-	float xb = 0.0f;
-	float xc = 0.0f;
-	float rog = rapp / ti.tgasm;
-	/* The if below was added to avoid ground reflectances lower than
-	 zero when ainr(1,1) greater than rapp/tgasm
-	 In such case either the choice of atmospheric model was not
-	 adequate for that image or the calculated apparent reflectance
-	 was too low. Run the model again for other conditions.
-	 The lines below just decrease ainr(1,1)/tgasm to avoid too
-	 bright pixels in the image. Check the output file to see if that
-	 has happened. */
-
-	float decrfact = 1.0f;
-	if (rog < (ainrpix / ti.tgasm))
-	{
-		do
-		{
-			decrfact = decrfact - 0.1f;
-			ainrpix = decrfact * ainrpix;
-		}
-		while(rog < (ainrpix / ti.tgasm));
-	}
-
-	rog = (rog - ainrpix / ti.tgasm) / ti.sutott / ti.sdtott;
-	rog = rog / (1.f + rog * ti.sast);
-	xa = (float)M_PI * ti.sb / ti.xmus / ti.seb / ti.tgasm / ti.sutott / ti.sdtott;
-	xb = ti.srotot / ti.sutott / ti.sdtott / ti.tgasm;
-	xc = ti.sast;
-
-	if (rog > 1) rog = 1;
-	if (rog < 0) rog = 0;
+    float xa = 0.0f;
+    float xb = 0.0f;
+    float xc = 0.0f;
+    float rog = rapp / ti.tgasm;
+    /* The if below was added to avoid ground reflectances lower than
+       zero when ainr(1,1) greater than rapp/tgasm
+       In such case either the choice of atmospheric model was not
+       adequate for that image or the calculated apparent reflectance
+       was too low. Run the model again for other conditions.
+       The lines below just decrease ainr(1,1)/tgasm to avoid too
+       bright pixels in the image. Check the output file to see if that
+       has happened. */
+
+    float decrfact = 1.0f;
+    if (rog < (ainrpix / ti.tgasm))
+    {
+	do
+	{
+	    decrfact = decrfact - 0.1f;
+	    ainrpix = decrfact * ainrpix;
+	}
+	while(rog < (ainrpix / ti.tgasm));
+    }
+
+    rog = (rog - ainrpix / ti.tgasm) / ti.sutott / ti.sdtott;
+    rog = rog / (1.f + rog * ti.sast);
+    xa = (float)M_PI * ti.sb / ti.xmus / ti.seb / ti.tgasm / ti.sutott / ti.sdtott;
+    xb = ti.srotot / ti.sutott / ti.sdtott / ti.tgasm;
+    xc = ti.sast;
+
+    if (rog > 1) rog = 1;
+    if (rog < 0) rog = 0;
 
     return rog;
 }

imagery/i.atcorr/gauss.cpp

@@ -10,17 +10,17 @@ float Gauss::angphi[13] = { 0.f, 30.f, 60.f, 90.f, 120.f, 150.f, 180.f, 210.f, 2
 /*  preliminary computations for gauss integration */
 void Gauss::init()
 {
-	int j;   
+    int j;   
 
-	/* convert angphi and angmu to radians */
-	for (j = 0; j < 13; ++j) angphi[j] = (float)(angphi[j] * M_PI / 180.f);
-	for (j = 0; j < 10; ++j) angmu[j] =	 (float)cos(angmu[j] * M_PI / 180.f);
+    /* convert angphi and angmu to radians */
+    for (j = 0; j < 13; ++j) angphi[j] = (float)(angphi[j] * M_PI / 180.f);
+    for (j = 0; j < 10; ++j) angmu[j] =	 (float)cos(angmu[j] * M_PI / 180.f);
 
-	/* calculate rm & gb */
+    /* calculate rm & gb */
 	
-	float anglem[mu2];
-	float weightm[mu2];
-	gauss (-1.f, 1.f, anglem, weightm, mu2);
+    float anglem[mu2];
+    float weightm[mu2];
+    gauss (-1.f, 1.f, anglem, weightm, mu2);
 
     gb[STDI(-mu)]   = 0;
     gb[STDI(0)]     = 0;
@@ -28,60 +28,60 @@ void Gauss::init()
     rm[STDI(-mu)]   = 0;
     rm[STDI(0)]     = 0;
     rm[STDI(mu)]    = 0;
-	/* do shift into rm & gb */
-	for (j = -mu+1; j <= -1; ++j)
+    /* do shift into rm & gb */
+    for (j = -mu+1; j <= -1; ++j)
     {
-      rm[-j] = anglem[mu + j - 1];
-      gb[-j] = weightm[mu + j - 1];
+	rm[-j] = anglem[mu + j - 1];
+	gb[-j] = weightm[mu + j - 1];
     }
 
-	for (j = 1; j <= mu-1; ++j)
+    for (j = 1; j <= mu-1; ++j)
     {
-      rm[2*mu - j] = anglem[mu + j - 2];
-      gb[2*mu - j] = weightm[mu + j - 2];
+	rm[2*mu - j] = anglem[mu + j - 2];
+	gb[2*mu - j] = weightm[mu + j - 2];
     }
 
-	/* calculate rp & gp */
-	gauss (0.f, (float)2 * M_PI, rp, gp, np);
+    /* calculate rp & gp */
+    gauss (0.f, (float)2 * M_PI, rp, gp, np);
 }
 
 
 /*	Compute for a given n, the gaussian quadrature (the n gaussian angles and the
-	their respective weights). The gaussian quadrature is used in numerical integration involving the
-	cosine of emergent or incident direction zenith angle. */
+  their respective weights). The gaussian quadrature is used in numerical integration involving the
+  cosine of emergent or incident direction zenith angle. */
 void Gauss::gauss (float a, float b, float *x, float *w, long int n)
 {
-	int m = (n + 1) / 2;
-	double xm = (b + a) / 2;
-	double xl = (b - a) / 2;
+    int m = (n + 1) / 2;
+    double xm = (b + a) / 2;
+    double xl = (b - a) / 2;
 
-	for(int i = 0; i < m; i++)
-	{
-		double 
-			z1, 
-			pp, 
-			z = cos(M_PI * (i + 0.75) / (n + 0.5));
+    for(int i = 0; i < m; i++)
+    {
+	double 
+	    z1, 
+	    pp, 
+	    z = cos(M_PI * (i + 0.75) / (n + 0.5));
 
-		do {
-			double p1 = 1;
-			double p2 = 0;
+	do {
+	    double p1 = 1;
+	    double p2 = 0;
 
-			for(int j = 0; j < n; j++)
-			{
-				double p3 = p2;
-				p2 = p1;
-				p1 = ((2 * j + 1) * z * p2 - j * p3) / (j+1);
-			}
+	    for(int j = 0; j < n; j++)
+	    {
+		double p3 = p2;
+		p2 = p1;
+		p1 = ((2 * j + 1) * z * p2 - j * p3) / (j+1);
+	    }
 
-			pp = n * (z * p1 - p2) / (z * z - 1);
-			z1 = z;
-			z = z1 - p1 / pp;
-		} while(fabs(z - z1) > 3e-14);
+	    pp = n * (z * p1 - p2) / (z * z - 1);
+	    z1 = z;
+	    z = z1 - p1 / pp;
+	} while(fabs(z - z1) > 3e-14);
 
-		if (fabs(z) < 3e-14) z = 0;
-		x[i] = (float)(xm - xl * z);
+	if (fabs(z) < 3e-14) z = 0;
+	x[i] = (float)(xm - xl * z);
         x[n - 1 - i] = (float)(xm + xl * z);
-		w[i] = (float)(2 * xl / ((1 - z * z) * pp * pp));
-		w[n - 1 - i] = w[i];
-	}
+	w[i] = (float)(2 * xl / ((1 - z * z) * pp * pp));
+	w[n - 1 - i] = w[i];
+    }
 }

imagery/i.atcorr/main.cpp

@@ -22,7 +22,7 @@
     supplying an elevation map has not been run to completion, because it
     takes to long and no sensible data for the test data was at hand.
     Testing would be welcomed. :)  
- ***************************************************************************/
+***************************************************************************/
 
 #include <stdlib.h>
 #include <math.h>
@@ -41,19 +41,19 @@ extern "C" {
 struct Options
 {
     /* options */
-	struct Option *iimg;    /* input satelite image */
-	struct Option *iscl;    /* input data is scaled to this range */
-	struct Option *ialt;    /* an input elevation map in km used to increase */
+    struct Option *iimg;    /* input satelite image */
+    struct Option *iscl;    /* input data is scaled to this range */
+    struct Option *ialt;    /* an input elevation map in km used to increase */
                             /* atmospheric correction accuracy, including this */
                             /* will make computations take much, much longer */
     struct Option *ivis;    /* an input visibility map in km (same purpose and effect as ialt) */
-	struct Option *icnd;    /* the input conditions file */
-	struct Option *oimg;    /* output image name */
-	struct Option *oscl;    /* scale the output data (reflectance values) to this range */
+    struct Option *icnd;    /* the input conditions file */
+    struct Option *oimg;    /* output image name */
+    struct Option *oscl;    /* scale the output data (reflectance values) to this range */
 
     /* flags */
-	struct Flag *oflt;      /* output data as floating point and do not round */
-	struct Flag *irad;      /* treat input values as reflectance instead of radiance values */
+    struct Flag *oflt;      /* output data as floating point and do not round */
+    struct Flag *irad;      /* treat input values as reflectance instead of radiance values */
     struct Flag *etmafter;  /* treat input data as a satelite image of type etm+ taken after July 1, 2000 */
     struct Flag *etmbefore; /* treat input data as a satelite image of type etm+ taken before July 1, 2000 */
     struct Flag *optimize;
@@ -81,42 +81,42 @@ static void read_scale (Option *, ScaleRange *);
 
 
 /* 
- Adjust the region to that of the input raster.
- Atmospheric corrections should be done on the whole
- satelite image, not just portions.
- */
+   Adjust the region to that of the input raster.
+   Atmospheric corrections should be done on the whole
+   satelite image, not just portions.
+*/
 static void adjust_region (char *name, char *mapset)
 {
-	struct Cell_head iimg_head;	/* the input image header file */
+    struct Cell_head iimg_head;	/* the input image header file */
 
-	if(G_get_cellhd(name, mapset, &iimg_head) < 0) 
-		G_fatal_error ("Unable to retreive header dat for input image");
+    if(G_get_cellhd(name, mapset, &iimg_head) < 0) 
+	G_fatal_error ("Unable to retreive header dat for input image");
 
-	if(G_set_window(&iimg_head) < 0) 
-		G_fatal_error ("Invalid graphics region coordinates");
+    if(G_set_window(&iimg_head) < 0) 
+	G_fatal_error ("Invalid graphics region coordinates");
 }
 
 
 /* Rounds a floating point cell value */
 static CELL round_c (FCELL x)
 {
-	if (x >= 0.0)
-		return (CELL)(x + .5);
+    if (x >= 0.0)
+	return (CELL)(x + .5);
 
-	return (CELL)(-(-x + .5));
+    return (CELL)(-(-x + .5));
 }
 
 
 /* Converts the buffer to cell and write it to disk */
 static void write_fp_to_cell (int ofd, FCELL* buf)
 {
-	CELL* cbuf;
-	int col;
+    CELL* cbuf;
+    int col;
 
-	cbuf = (CELL*)G_allocate_raster_buf(CELL_TYPE);
+    cbuf = (CELL*)G_allocate_raster_buf(CELL_TYPE);
 
-	for(col = 0; col < G_window_cols(); col++) cbuf[col] = round_c(buf[col]);
-	G_put_raster_row(ofd, cbuf, CELL_TYPE);
+    for(col = 0; col < G_window_cols(); col++) cbuf[col] = round_c(buf[col]);
+    G_put_raster_row(ofd, cbuf, CELL_TYPE);
 }
 
 
@@ -133,85 +133,85 @@ class TICache
 public:
     TICache() { p = 0; for(int i = 0; i < MAX_TIs; i++) alts[i] = -1; }
     int search(float alt) { 
-	    for(int i = 0; i < MAX_TIs; i++) 
-		    if(alt == alts[i]) 
-		    {
-			    hit++;
-			    return i;
-		    } 
-	    mis++;
-	    return -1; 
+	for(int i = 0; i < MAX_TIs; i++) 
+	    if(alt == alts[i]) 
+	    {
+		hit++;
+		return i;
+	    } 
+	mis++;
+	return -1; 
     }
 
     TransformInput get(int i) { return tis[i]; }
     void add(TransformInput ti, float alt) { 
-	    tis[p] = ti; 
-	    alts[p] = alt; 
-	    p++; 
-	    if(p >= MAX_TIs) p = 0; 
+	tis[p] = ti; 
+	alts[p] = alt; 
+	p++; 
+	if(p >= MAX_TIs) p = 0; 
     }
 };
 
 
 /* the transform input map, is a array of ticaches.
- The first key is the visibility which matches to a TICache for the altitudes.
- This code is horrible, i just spent 20min writing and 5min debugging it. */
+   The first key is the visibility which matches to a TICache for the altitudes.
+   This code is horrible, i just spent 20min writing and 5min debugging it. */
 class TIMap
 {
-	enum TIMapSize
-	{
-		MAX_TICs = 128 /* this value is a guess. It means that 1024 TI's will be the max combinations of vis/alt pairs */
-	};
+    enum TIMapSize
+    {
+	MAX_TICs = 128 /* this value is a guess. It means that 1024 TI's will be the max combinations of vis/alt pairs */
+    };
 
-	TICache tic[MAX_TICs]; /* array of TICaches */
-	float visi[MAX_TICs];
-	int p;
+    TICache tic[MAX_TICs]; /* array of TICaches */
+    float visi[MAX_TICs];
+    int p;
 
 public:
-	struct Position
-	{
-		int i, j;
-		Position() : i(-1), j(-1) {}
-		Position(int x, int y) : i(x), j(y) {}
-		bool valid() { return i != -1 && j != -1; }
-	};	
+    struct Position
+    {
+	int i, j;
+	Position() : i(-1), j(-1) {}
+	Position(int x, int y) : i(x), j(y) {}
+	bool valid() { return i != -1 && j != -1; }
+    };	
 	
-	TIMap() { p = 0; for(int i = 0; i < MAX_TICs; i++) visi[i] = -1; }
-	Position search(float vis, float alt) { 
-		for(int i = 0; i < MAX_TICs; i++)
-			if(vis == visi[i]) {
-				Position pos;
-				pos.i = i;
-				pos.j = tic[i].search(alt);
-				return pos;
-			} 
-		return Position();
-	}
+    TIMap() { p = 0; for(int i = 0; i < MAX_TICs; i++) visi[i] = -1; }
+    Position search(float vis, float alt) { 
+	for(int i = 0; i < MAX_TICs; i++)
+	    if(vis == visi[i]) {
+		Position pos;
+		pos.i = i;
+		pos.j = tic[i].search(alt);
+		return pos;
+	    } 
+	return Position();
+    }
 
-	TransformInput get(Position pos) { return tic[pos.i].get(pos.j); }
+    TransformInput get(Position pos) { return tic[pos.i].get(pos.j); }
 
-	void add(TransformInput ti, float vis, float alt) {
-		tic[p].add(ti, alt);
-		visi[p] = vis;
-		p++;
-		if(p >= MAX_TICs) p = 0;
-	}
+    void add(TransformInput ti, float vis, float alt) {
+	tic[p].add(ti, alt);
+	visi[p] = vis;
+	p++;
+	if(p >= MAX_TICs) p = 0;
+    }
 };
 
 
 struct IntPair
 {
-	FCELL x;
-	FCELL y;
+    FCELL x;
+    FCELL y;
 	
-	IntPair(FCELL i, FCELL j) : x(i), y(j) {}
+    IntPair(FCELL i, FCELL j) : x(i), y(j) {}
 	
-	bool operator<(const IntPair& b) const
+    bool operator<(const IntPair& b) const
 	{
-		if(x < b.x) return true;
-		else if(x > b.x) return false;
-		else if(y < b.y) return true;
-		return false;
+	    if(x < b.x) return true;
+	    else if(x > b.x) return false;
+	    else if(y < b.y) return true;
+	    return false;
 	}	
 };
 
@@ -221,51 +221,51 @@ typedef std::map<IntPair, TransformInput> CacheMap;
 
 const TransformInput& optimize_va (const FCELL& vis, const FCELL& alt)
 {
-	static CacheMap timap;
-	static TransformInput ti;
+    static CacheMap timap;
+    static TransformInput ti;
 
-	IntPair key(vis, alt);
-	CacheMap::iterator it = timap.find(key);
+    IntPair key(vis, alt);
+    CacheMap::iterator it = timap.find(key);
 
-	if(it != timap.end()) /* search found key */
-	{
-		ti = (*it).second;
-	}
-	else
-	{
-		pre_compute_hv(alt, vis);
-		ti = compute();
-		timap.insert(std::make_pair(key, ti));
-	}
+    if(it != timap.end()) /* search found key */
+    {
+	ti = (*it).second;
+    }
+    else
+    {
+	pre_compute_hv(alt, vis);
+	ti = compute();
+	timap.insert(std::make_pair(key, ti));
+    }
 	
-	return ti;
+    return ti;
 }	
 
 
 /* Process the raster and do atmospheric corrections.
-Params:
- * INPUT FILE
- ifd: input file descriptor
- iref: input file has radiance values (default is reflectance) ?
- iscale: input file's range (default is min = 0, max = 255)
- ialt_fd: height map file descriptor, negative if global value is used
- ivis_fd: visibility map file descriptor, negative if global value is used
-
- * OUTPUT FILE
- ofd: output file descriptor
- oflt: if true use FCELL_TYPE for output
- oscale: output file's range (default is min = 0, max = 255)
+   Params:
+   * INPUT FILE
+   ifd: input file descriptor
+   iref: input file has radiance values (default is reflectance) ?
+   iscale: input file's range (default is min = 0, max = 255)
+   ialt_fd: height map file descriptor, negative if global value is used
+   ivis_fd: visibility map file descriptor, negative if global value is used
+
+   * OUTPUT FILE
+   ofd: output file descriptor
+   oflt: if true use FCELL_TYPE for output
+   oscale: output file's range (default is min = 0, max = 255)
 */
 static void process_raster (int ifd, InputMask imask, ScaleRange iscale,
-                int ialt_fd, int ivis_fd, int ofd, bool oflt,
-                ScaleRange oscale, bool optimize)
+			    int ialt_fd, int ivis_fd, int ofd, bool oflt,
+			    ScaleRange oscale, bool optimize)
 {
-	FCELL* buf;         /* buffer for the input values */
+    FCELL* buf;         /* buffer for the input values */
     FCELL* alt = NULL;         /* buffer for the elevation values */
     FCELL* vis = NULL;         /* buffer for the visibility values */
     FCELL  prev_alt = -1.f;
     FCELL  prev_vis = -1.f;
-	int row, col;
+    int row, col;
 
     /* do initial computation with global elevation and visibility values */
     TransformInput ti;
@@ -274,36 +274,36 @@ static void process_raster (int ifd, InputMask imask, ScaleRange iscale,
     TICache ticache;    /* use this to increase computation speed when an elevation map with categories are given */
 	
     /* allocate memory for buffers */
-	buf = (FCELL*)G_allocate_raster_buf(FCELL_TYPE);
+    buf = (FCELL*)G_allocate_raster_buf(FCELL_TYPE);
     if(ialt_fd >= 0) alt = (FCELL*)G_allocate_raster_buf(FCELL_TYPE);
     if(ivis_fd >= 0) vis = (FCELL*)G_allocate_raster_buf(FCELL_TYPE);
 
     fprintf(stderr, "Percent complete: ");
 
-	for(row = 0; row < G_window_rows(); row++)
-	{
-        	G_percent(row, G_window_rows(), 1);     /* keep the user informed of our progress */
+    for(row = 0; row < G_window_rows(); row++)
+    {
+	G_percent(row, G_window_rows(), 1);     /* keep the user informed of our progress */
 		
         /* read the next row */
-		if(G_get_raster_row(ifd, buf, row, FCELL_TYPE) < 0)
-			G_fatal_error ("Unable to read from input file");
+	if(G_get_raster_row(ifd, buf, row, FCELL_TYPE) < 0)
+	    G_fatal_error ("Unable to read from input file");
 
         /* read the next row of elevation values */
         if(ialt_fd >= 0)
-    		if(G_get_raster_row(ialt_fd, alt, row, FCELL_TYPE) < 0)
-	    		G_fatal_error ("Unable to read from elevation raster");
+	    if(G_get_raster_row(ialt_fd, alt, row, FCELL_TYPE) < 0)
+		G_fatal_error ("Unable to read from elevation raster");
 
         /* read the next row of elevation values */
         if(ivis_fd >= 0)
-    		if(G_get_raster_row(ivis_fd, vis, row, FCELL_TYPE) < 0)
-	    		G_fatal_error ("Unable to read from visibility raster");
+	    if(G_get_raster_row(ivis_fd, vis, row, FCELL_TYPE) < 0)
+		G_fatal_error ("Unable to read from visibility raster");
 
         /* loop over all the values in the row */
-		for(col = 0; col < G_window_cols(); col++)
-		{
+	for(col = 0; col < G_window_cols(); col++)
+	{
 /* TODO: use G_set_f_null_value()?? */
-		if(vis && isnan(vis[col]) || alt && isnan(alt[col]) || isnan(buf[col])) {buf[col] = FP_NAN; continue;}
-      		alt[col] /= 1000.0f; /* converting to km from input which should be in meter */
+	    if(vis && isnan(vis[col]) || alt && isnan(alt[col]) || isnan(buf[col])) {buf[col] = FP_NAN; continue;}
+	    alt[col] /= 1000.0f; /* converting to km from input which should be in meter */
 
             /* check if both maps are active and if whether any value has changed */
             if((ialt_fd >= 0) && (ivis_fd >= 0) && ((prev_vis != vis[col]) || (prev_alt != alt[col])))
@@ -312,8 +312,8 @@ static void process_raster (int ifd, InputMask imask, ScaleRange iscale,
                	prev_vis = vis[col];
  		if(optimize) ti = optimize_va(vis[col], alt[col]); /* try to optimize? */
 		else { /* no optimizations */
-		   pre_compute_hv(alt[col], vis[col]);
-               	   ti = compute();
+		    pre_compute_hv(alt[col], vis[col]);
+		    ti = compute();
 		}	
             }
             else    /* only one of the maps is being used */
@@ -372,16 +372,16 @@ static void process_raster (int ifd, InputMask imask, ScaleRange iscale,
             buf[col] = buf[col] * ((float)oscale.max - (float)oscale.min) + oscale.min;
 
             if(~oflt && (buf[col] > (float)oscale.max))
-              G_warning ("The output data will overflow. Reflectance > 100%%");
-		}
+		G_warning ("The output data will overflow. Reflectance > 100%%");
+	}
 
         /* write output */
-		if(oflt) G_put_raster_row(ofd, buf, FCELL_TYPE);
-		else write_fp_to_cell(ofd, buf);
-	}
+	if(oflt) G_put_raster_row(ofd, buf, FCELL_TYPE);
+	else write_fp_to_cell(ofd, buf);
+    }
 
     /* free allocated memory */
-	G_free(buf);
+    G_free(buf);
     if(ialt_fd >= 0) G_free(alt);
     if(ivis_fd >= 0) G_free(vis);
 }
@@ -391,111 +391,111 @@ static void process_raster (int ifd, InputMask imask, ScaleRange iscale,
 /* Copy the colors from map named iname to the map named oname */
 static void copy_colors (char *iname, char *imapset, char *oname)
 {
-	struct Colors colors;
+    struct Colors colors;
 
-	G_read_colors(iname, imapset, &colors);
-	G_write_colors(oname, G_mapset(), &colors);
+    G_read_colors(iname, imapset, &colors);
+    G_write_colors(oname, G_mapset(), &colors);
 }
 
 
 /* Define our module so that Grass can print it if the user wants to know more. */
 static void define_module (void)
 {
-	struct GModule *module;
-
-	module = G_define_module();
-	module->label = _("Performs atmospheric correction using the 6S algorithm.");
-	module->description =
-	 _("6S - Second Simulation of Satellite Signal in the Solar Spectrum.");
-	/* 
-	 " Incorporated into Grass by Christo A. Zietsman, January 2003.\n"
-	 " Converted from Fortran to C by Christo A. Zietsman, November 2002.\n\n"
-	 " Adapted by Mauro A. Homem Antunes for atmopheric corrections of\n"
-	 " remotely sensed images in raw format (.RAW) of 8 bits.\n"
-	 " April 4, 2001.\n\n"
-	 " Please refer to the following paper and acknowledge the authors of\n"
-	 " the model:\n"
-	 " Vermote, E.F., Tanre, D., Deuze, J.L., Herman, M., and Morcrette,\n"
-	 "    J.J., (1997), Second simulation of the satellite signal in\n"
-	 "    the solar spectrum, 6S: An overview., IEEE Trans. Geosc.\n"
-	 "    and Remote Sens. 35(3):675-686.\n"
-	 " The code is provided as is and is not to be sold. See notes on\n"
-	 " http://loasys.univ-lille1.fr/informatique/sixs_gb.html\n"
-	 " http://www.ltid.inpe.br/dsr/mauro/6s/index.html\n"
-	 " and on http://www.cs.sun.ac.za/~caz/index.html\n";*/
+    struct GModule *module;
+
+    module = G_define_module();
+    module->label = _("Performs atmospheric correction using the 6S algorithm.");
+    module->description =
+	_("6S - Second Simulation of Satellite Signal in the Solar Spectrum.");
+    /* 
+       " Incorporated into Grass by Christo A. Zietsman, January 2003.\n"
+       " Converted from Fortran to C by Christo A. Zietsman, November 2002.\n\n"
+       " Adapted by Mauro A. Homem Antunes for atmopheric corrections of\n"
+       " remotely sensed images in raw format (.RAW) of 8 bits.\n"
+       " April 4, 2001.\n\n"
+       " Please refer to the following paper and acknowledge the authors of\n"
+       " the model:\n"
+       " Vermote, E.F., Tanre, D., Deuze, J.L., Herman, M., and Morcrette,\n"
+       "    J.J., (1997), Second simulation of the satellite signal in\n"
+       "    the solar spectrum, 6S: An overview., IEEE Trans. Geosc.\n"
+       "    and Remote Sens. 35(3):675-686.\n"
+       " The code is provided as is and is not to be sold. See notes on\n"
+       " http://loasys.univ-lille1.fr/informatique/sixs_gb.html\n"
+       " http://www.ltid.inpe.br/dsr/mauro/6s/index.html\n"
+       " and on http://www.cs.sun.ac.za/~caz/index.html\n";*/
 }
 
 
 /* Define the options and flags */
 static struct Options define_options (void)
 {
-	struct Options opts;
+    struct Options opts;
 
-	opts.iimg = G_define_standard_option (G_OPT_R_INPUT);
-	opts.iimg->key		= "iimg";
-	opts.iimg->description	= "Input imagery map to be corrected";
+    opts.iimg = G_define_standard_option (G_OPT_R_INPUT);
+    opts.iimg->key		= "iimg";
+    opts.iimg->description	= "Input imagery map to be corrected";
 /*	opts.iimg->answer	= "ETM4_400x400.raw"; */
 
-	opts.iscl = G_define_option();
-	opts.iscl->key          = "iscl";
-	opts.iscl->type         = TYPE_INTEGER;
-	opts.iscl->key_desc     = "Input scale range";
-	opts.iscl->required     = NO;
-	opts.iscl->answer       = "0,255";
-	opts.iscl->description  = "Input imagery range [0,255]";
-
-	opts.ialt = G_define_standard_option (G_OPT_R_INPUT);
-	opts.ialt->key		= "ialt";
-	opts.ialt->required	= NO;
-	opts.ialt->answer	= "dem_float";
-	opts.ialt->description	= "Input altitude map in m (optional)";
-
-	opts.ivis = G_define_standard_option (G_OPT_R_INPUT);
-	opts.ivis->key		= "ivis";
-	opts.ivis->required	= NO;
+    opts.iscl = G_define_option();
+    opts.iscl->key          = "iscl";
+    opts.iscl->type         = TYPE_INTEGER;
+    opts.iscl->key_desc     = "Input scale range";
+    opts.iscl->required     = NO;
+    opts.iscl->answer       = "0,255";
+    opts.iscl->description  = "Input imagery range [0,255]";
+
+    opts.ialt = G_define_standard_option (G_OPT_R_INPUT);
+    opts.ialt->key		= "ialt";
+    opts.ialt->required	= NO;
+    opts.ialt->answer	= "dem_float";
+    opts.ialt->description	= "Input altitude map in m (optional)";
+
+    opts.ivis = G_define_standard_option (G_OPT_R_INPUT);
+    opts.ivis->key		= "ivis";
+    opts.ivis->required	= NO;
 /*	opts.ivis->answer	= "visibility"; */
-	opts.ivis->description	= "Input visibility map in km (optional)";
+    opts.ivis->description	= "Input visibility map in km (optional)";
 
-	opts.icnd = G_define_standard_option (G_OPT_F_INPUT);
-	opts.icnd->key		= "icnd";
-	opts.icnd->required	= YES;
+    opts.icnd = G_define_standard_option (G_OPT_F_INPUT);
+    opts.icnd->key		= "icnd";
+    opts.icnd->required	= YES;
 /*	opts.icnd->answer	= "ETM4_atmospheric_input_GRASS.txt"; */
-	opts.icnd->description	= "6S input text file";
+    opts.icnd->description	= "6S input text file";
 
-	opts.oimg = G_define_standard_option (G_OPT_R_OUTPUT);
-	opts.oimg->key		= "oimg";
+    opts.oimg = G_define_standard_option (G_OPT_R_OUTPUT);
+    opts.oimg->key		= "oimg";
 /*	opts.oimg->answer	= "6s_output_file"; */
-	opts.oimg->description	= "6S output imagery map";
+    opts.oimg->description	= "6S output imagery map";
 
-	opts.oscl = G_define_option();
-	opts.oscl->key          = "oscl";
-	opts.oscl->type         = TYPE_INTEGER;
-	opts.oscl->key_desc     = "Output scale range";
-	opts.oscl->required     = YES;
-	opts.oscl->answer       = "0,255";
-	opts.oscl->description  = "Rescale output imagery map [0,255]";
+    opts.oscl = G_define_option();
+    opts.oscl->key          = "oscl";
+    opts.oscl->type         = TYPE_INTEGER;
+    opts.oscl->key_desc     = "Output scale range";
+    opts.oscl->required     = YES;
+    opts.oscl->answer       = "0,255";
+    opts.oscl->description  = "Rescale output imagery map [0,255]";
 
-	opts.oflt = G_define_flag();
-	opts.oflt->key = 'f';
-	opts.oflt->description = "Output raster is floating point";
+    opts.oflt = G_define_flag();
+    opts.oflt->key = 'f';
+    opts.oflt->description = "Output raster is floating point";
 
-	opts.irad = G_define_flag();
-	opts.irad->key = 'r';
-	opts.irad->description = "Input map converted to reflectance (default is radiance)";
+    opts.irad = G_define_flag();
+    opts.irad->key = 'r';
+    opts.irad->description = "Input map converted to reflectance (default is radiance)";
 
-	opts.etmafter = G_define_flag();
-	opts.etmafter->key = 'a';
-	opts.etmafter->description = "Input from ETM+ image taken after July 1, 2000";
+    opts.etmafter = G_define_flag();
+    opts.etmafter->key = 'a';
+    opts.etmafter->description = "Input from ETM+ image taken after July 1, 2000";
 
-	opts.etmbefore = G_define_flag();
-	opts.etmbefore->key = 'b';
-	opts.etmbefore->description = "Input from ETM+ image taken before July 1, 2000";
+    opts.etmbefore = G_define_flag();
+    opts.etmbefore->key = 'b';
+    opts.etmbefore->description = "Input from ETM+ image taken before July 1, 2000";
 
-	opts.optimize = G_define_flag();
-	opts.optimize->key = 'o';
-	opts.optimize->description = "Try to increase computation speed when categorized altitude or/and visibility map is used.";
+    opts.optimize = G_define_flag();
+    opts.optimize->key = 'o';
+    opts.optimize->description = "Try to increase computation speed when categorized altitude or/and visibility map is used.";
 
-	return opts;
+    return opts;
 }
 
 /* Read the min and max values from the iscl and oscl options */
@@ -532,88 +532,88 @@ static void read_scale (Option *scl, ScaleRange &range)
 
 int main(int argc, char* argv[])
 {
-	struct Options opts;        
+    struct Options opts;        
     struct ScaleRange iscale;   /* input file's data is scaled to this interval */
     struct ScaleRange oscale;   /* output file's scale */
-	int iimg_fd;	        /* input image's file descriptor */
-	int oimg_fd;	        /* output image's file descriptor */
-	int ialt_fd = -1;       /* input elevation map's file descriptor */
+    int iimg_fd;	        /* input image's file descriptor */
+    int oimg_fd;	        /* output image's file descriptor */
+    int ialt_fd = -1;       /* input elevation map's file descriptor */
     int ivis_fd = -1;       /* input visibility map's file descriptor */
     char *iimg_mapset, *ialt_mapset, *iviz_mapset;
     
-	/* Define module */
-	define_module();
+    /* Define module */
+    define_module();
   
     /* Define the different input options */
-	opts = define_options();
+    opts = define_options();
 
-	/**** Start ****/
-	G_gisinit(argv[0]);
-	if (G_parser(argc, argv) < 0)
-		exit (EXIT_FAILURE);
+    /**** Start ****/
+    G_gisinit(argv[0]);
+    if (G_parser(argc, argv) < 0)
+	exit (EXIT_FAILURE);
 
-	/* open input raster */
-	if ( (iimg_mapset = G_find_cell2 ( opts.iimg->answer, "") ) == NULL )
-	     G_fatal_error ( _("Raster map <%s> not found"), opts.iimg->answer);
-	if((iimg_fd = G_open_cell_old(opts.iimg->answer, iimg_mapset)) < 0)
-		G_fatal_error ("Unable to open input raster");
+    /* open input raster */
+    if ( (iimg_mapset = G_find_cell2 ( opts.iimg->answer, "") ) == NULL )
+	G_fatal_error ( _("Raster map <%s> not found"), opts.iimg->answer);
+    if((iimg_fd = G_open_cell_old(opts.iimg->answer, iimg_mapset)) < 0)
+	G_fatal_error ("Unable to open input raster");
 
-	adjust_region(opts.iimg->answer, iimg_mapset);
+    adjust_region(opts.iimg->answer, iimg_mapset);
         
-        if(opts.ialt->answer) {
-	  if ( (ialt_mapset = G_find_cell2 ( opts.ialt->answer, "") ) == NULL )
+    if(opts.ialt->answer) {
+	if ( (ialt_mapset = G_find_cell2 ( opts.ialt->answer, "") ) == NULL )
 	    G_fatal_error ( _("Raster map <%s> not found"), opts.ialt->answer);
-	  if((ialt_fd = G_open_cell_old(opts.ialt->answer, ialt_mapset)) < 0)
+	if((ialt_fd = G_open_cell_old(opts.ialt->answer, ialt_mapset)) < 0)
             G_warning ("Unable to open DEM raster");
-	}
+    }
 
-	if(opts.ivis->answer) {
-	  if ( (iviz_mapset = G_find_cell2 ( opts.ivis->answer, "") ) == NULL )
-	       G_fatal_error ( _("Raster map <%s> not found"), opts.ivis->answer);
-          if((ivis_fd = G_open_cell_old(opts.ivis->answer, iviz_mapset)) < 0)
+    if(opts.ivis->answer) {
+	if ( (iviz_mapset = G_find_cell2 ( opts.ivis->answer, "") ) == NULL )
+	    G_fatal_error ( _("Raster map <%s> not found"), opts.ivis->answer);
+	if((ivis_fd = G_open_cell_old(opts.ivis->answer, iviz_mapset)) < 0)
             G_warning ("Unable to open visibility raster");
-	}
+    }
                 
-	/* open a floating point raster or not? */
-	if(opts.oflt->answer)
-	{
-		if((oimg_fd = G_open_fp_cell_new(opts.oimg->answer)) < 0)
-			G_fatal_error ("Unable to create output raster");
-	}
-	else
-	{
-		if((oimg_fd = G_open_raster_new(opts.oimg->answer, CELL_TYPE)) < 0)
-			G_fatal_error ("Unable to create output raster");
-	}
+    /* open a floating point raster or not? */
+    if(opts.oflt->answer)
+    {
+	if((oimg_fd = G_open_fp_cell_new(opts.oimg->answer)) < 0)
+	    G_fatal_error ("Unable to create output raster");
+    }
+    else
+    {
+	if((oimg_fd = G_open_raster_new(opts.oimg->answer, CELL_TYPE)) < 0)
+	    G_fatal_error ("Unable to create output raster");
+    }
 
     /* read the scale parameters */
     read_scale(opts.iscl, iscale);
     read_scale(opts.oscl, oscale);
 
     /* initialize this 6s computation and parse the input conditions file */
-	init_6S(opts.icnd->answer);
+    init_6S(opts.icnd->answer);
 	
     InputMask imask = RADIANCE;         /* the input mask tells us what transformations if any
-                                         needs to be done to make our input values, reflectance
-                                         values scaled between 0 and 1 */
+					   needs to be done to make our input values, reflectance
+					   values scaled between 0 and 1 */
     if(opts.irad->answer) imask = REFLECTANCE;
     if(opts.etmbefore->answer) imask = (InputMask)(imask | ETM_BEFORE);
     if(opts.etmafter->answer) imask = (InputMask)(imask | ETM_AFTER);
 
     /* process the input raster and produce our atmospheric corrected output raster. */
-	process_raster(iimg_fd, imask, iscale, ialt_fd, ivis_fd,
+    process_raster(iimg_fd, imask, iscale, ialt_fd, ivis_fd,
                    oimg_fd, opts.oflt->answer, oscale, opts.optimize->answer);
 
 
     /* Close the input and output file descriptors */
-	G_close_cell(iimg_fd);
+    G_close_cell(iimg_fd);
     if(opts.ialt->answer) G_close_cell(ialt_fd);
     if(opts.ivis->answer) G_close_cell(ivis_fd);
-	G_close_cell(oimg_fd);
+    G_close_cell(oimg_fd);
 
     /* Copy the colors of the input raster to the output raster.
        Scaling is ignored and color ranges might not be correct. */
-	copy_colors(opts.iimg->answer, iimg_mapset, opts.oimg->answer);
+    copy_colors(opts.iimg->answer, iimg_mapset, opts.oimg->answer);
 
-	exit (EXIT_SUCCESS);
+    exit (EXIT_SUCCESS);
 }