i.atcorr: avoid float to double to float conversion, causing numerical instability

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

imagery/i.atcorr/6s.cpp

@@ -22,10 +22,10 @@ extern "C" {
 extern void discom(const GeomCond &geom, const AtmosModel &atms,
                    const AerosolModel &aero, const AerosolConcentration &aerocon,
                    const Altitude &alt, const IWave &iwave);
-extern void specinterp(const float wl, float& tamoy, float& tamoyp, float& pizmoy, float& pizmoyp,
+extern void specinterp(const double wl, double& tamoy, double& tamoyp, double& pizmoy, double& 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);
+extern void enviro (const double difr, const double difa, const double r, const double palt,
+		    const double xmuv, double& fra, double& fae, double& fr);
 void printOutput(); // forward declare this function so that it can be used in init_6S
 
 
@@ -91,14 +91,14 @@ int init_6S(char* icnd_name)
 	c**********************************************************************/
 
     /* NOTE: wlmoy is not affected by a height and/or vis change */
-    float wlmoy;
+    double wlmoy;
     if(iwave.iwave != -1) wlmoy = iwave.equivwl();
     else wlmoy = iwave.wl;
 
     iwave.wlmoy = wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
-    float tamoy, tamoyp, pizmoy, pizmoyp;
+    double tamoy, tamoyp, pizmoy, pizmoyp;
     if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
 
     printOutput();
@@ -107,45 +107,45 @@ int init_6S(char* icnd_name)
 }
 
 /* Only update those objects that are affected by a height and vis change */
-void pre_compute_hv(const float height, const float vis)
+void pre_compute_hv(const double height, const double vis)
 {
     atms = original_atms;
     aerocon.set_visibility(vis, atms);
     alt.set_height(height);
     alt.init(atms, aerocon);
    
-    float wlmoy = iwave.wlmoy;
+    double wlmoy = iwave.wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
-    float tamoy, tamoyp, pizmoy, pizmoyp;
+    double 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 */
-void pre_compute_v(const float vis)
+void pre_compute_v(const double vis)
 {
     atms = original_atms;
     aerocon.set_visibility(vis, atms);
     alt.init(atms, aerocon);
 
-    float wlmoy = iwave.wlmoy;
+    double wlmoy = iwave.wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
-    float tamoy, tamoyp, pizmoy, pizmoyp;
+    double 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 */
-void pre_compute_h(const float height)
+void pre_compute_h(const double height)
 {
     atms = original_atms;
     alt.set_height(height);
     alt.init(atms, aerocon);
 
-    float wlmoy = iwave.wlmoy;
+    double wlmoy = iwave.wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
-    float tamoy, tamoyp, pizmoy, pizmoyp;
+    double tamoy, tamoyp, pizmoy, pizmoyp;
     if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
 }
 
@@ -194,7 +194,7 @@ void printOutput()
 	string(" spectral volcanic debris reflectance  ")
     };
 
-    float rocave = 0;       /* block of code in Fortran will always compute 0 */
+    double rocave = 0;       /* block of code in Fortran will always compute 0 */
     ostringstream s;
     s.setf(ios::fixed, ios::floatfield);
     s << setprecision(3);
@@ -241,71 +241,71 @@ void printOutput()
 
 TransformInput compute()
 {
-    const float accu3 = 1e-07;
+    const double 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];
+    double fr = 0;
+    double rad = 0;
+    double sb = 0;
+    double seb = 0;
+    double refet = 0;
+    double refet1 = 0;
+    double refet2 = 0;
+    double refet3 = 0;
+    double alumet = 0;
+    double tgasm = 0;
+    double rog = 0;
+    double dgasm = 0;
+    double ugasm = 0;
+    double sdwava = 0;
+    double sdozon = 0;
+    double sddica = 0;
+    double sdoxyg = 0;
+    double sdniox = 0;
+    double sdmoca = 0;
+    double sdmeth = 0;
+
+    double suwava = 0;
+    double suozon = 0;
+    double sudica = 0;
+    double suoxyg = 0;
+    double suniox = 0;
+    double sumoca = 0;
+    double sumeth = 0;
+    double stwava = 0;
+    double stozon = 0;
+    double stdica = 0;
+    double stoxyg = 0;
+    double stniox = 0;
+    double stmoca = 0;
+    double stmeth = 0;
+    double sodray = 0;
+    double sodrayp = 0;
+    double sodaer = 0;
+    double sodaerp = 0;
+    double sodtot = 0;
+    double sodtotp = 0;
+    double fophsr = 0;
+    double fophsa = 0;
+    double sroray = 0;
+    double sroaer = 0;
+    double srotot = 0;
+    double ssdaer = 0;
+    double sdtotr = 0;
+    double sdtota = 0;
+    double sdtott = 0;
+    double sutotr = 0;
+    double sutota = 0;
+    double sutott = 0;
+    double sasr = 0;
+    double sasa = 0;
+    double sast = 0;
+
+    double ani[2][3];
+    double aini[2][3];
+    double anr[2][3];
+    double ainr[2][3];
 
     for(i = 0; i < 2; i++)
 	for(j = 0; j < 3; j++)
@@ -327,24 +327,24 @@ TransformInput compute()
     int l;
     for(l = iwave.iinf; l <= iwave.isup; l++)
     {
-        float sbor = iwave.ffu.s[l];
+        double sbor = iwave.ffu.s[l];
 
         if(l == iwave.iinf || l == iwave.isup) sbor *= 0.5f;
         if(iwave.iwave == -1) sbor = 1.0f / step;
 
-        float roc = 0; /* rocl[l]; */
-        float roe = 0; /* roel[l]; */
-        float wl = 0.25f + l * step;
+        double roc = 0; /* rocl[l]; */
+        double roe = 0; /* roel[l]; */
+        double wl = 0.25f + l * step;
 
 	AbstraStruct as;
-	float uwus, uo3us;		/* initialized in abstra */
+	double uwus, uo3us;		/* initialized in abstra */
 
-	abstra(atms, alt, wl, (float)geom.xmus, (float)geom.xmuv, atms.uw / 2.0f, atms.uo3,
+	abstra(atms, alt, wl, (double)geom.xmus, (double)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;
+	double attwava = as.ttwava;
 
-	abstra(atms, alt, wl, (float)geom.xmus, (float)geom.xmuv, atms.uw, atms.uo3,
+	abstra(atms, alt, wl, (double)geom.xmus, (double)geom.xmuv, atms.uw, atms.uo3,
 	       uwus, uo3us, alt.puw, alt.puo3, alt.puwus, alt.puo3us, as);
 
         if (as.dtwava < accu3) as.dtwava = 0;
@@ -366,36 +366,36 @@ TransformInput compute()
         if (as.ttmeth < accu3) as.ttmeth = 0;
         if (as.ttmoca < accu3) as.ttmeth = 0;
 
-        float swl = iwave.solirr(wl);
+        double swl = iwave.solirr(wl);
         swl = swl * geom.dsol;
-        float coef = sbor * step * swl;
+        double 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);
+	interp(aero.iaer, alt.idatmp, wl, aerocon.taer55, alt.taer55p, (double)geom.xmud, is);
 
 
-        float dgtot = as.dtwava * as.dtozon * as.dtdica * as.dtoxyg * as.dtniox * as.dtmeth * as.dtmoca;
-        float tgtot = as.ttwava * as.ttozon * as.ttdica * as.ttoxyg * as.ttniox * as.ttmeth * as.ttmoca;
-        float ugtot = as.utwava * as.utozon * as.utdica * as.utoxyg * as.utniox * as.utmeth * as.utmoca;
-        float tgp1 = as.ttozon * as.ttdica * as.ttoxyg * as.ttniox * as.ttmeth * as.ttmoca;
-        float tgp2 = attwava * as.ttozon * as.ttdica * as.ttoxyg * as.ttniox * as.ttmeth * as.ttmoca;
-        float edifr = (float)(is.utotr - exp(-is.trayp / geom.xmuv));
-        float edifa = (float)(is.utota - exp(-is.taerp / geom.xmuv));
+        double dgtot = as.dtwava * as.dtozon * as.dtdica * as.dtoxyg * as.dtniox * as.dtmeth * as.dtmoca;
+        double tgtot = as.ttwava * as.ttozon * as.ttdica * as.ttoxyg * as.ttniox * as.ttmeth * as.ttmoca;
+        double ugtot = as.utwava * as.utozon * as.utdica * as.utoxyg * as.utniox * as.utmeth * as.utmoca;
+        double tgp1 = as.ttozon * as.ttdica * as.ttoxyg * as.ttniox * as.ttmeth * as.ttmoca;
+        double tgp2 = attwava * as.ttozon * as.ttdica * as.ttoxyg * as.ttniox * as.ttmeth * as.ttmoca;
+        double edifr = (double)(is.utotr - exp(-is.trayp / geom.xmuv));
+        double edifa = (double)(is.utota - exp(-is.taerp / geom.xmuv));
 
 
-	float fra, fae;
-	enviro(edifr, edifa, rad, alt.palt, (float)geom.xmuv, fra, fae, fr);
+	double fra, fae;
+	enviro(edifr, edifa, rad, alt.palt, (double)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)
+	double avr = roc * fr + (1 - fr) * roe;
+	double rsurf = (double)(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;
-        float romeas1 = ratm1 + rsurf * tgtot;
-        float romeas2 = ratm2 + rsurf * tgtot;
-        float romeas3 = ratm3 + rsurf * tgtot;
+        double ratm1 = (is.romix - is.rorayl) * tgtot + is.rorayl * tgp1;
+        double ratm3 = is.romix * tgp1;
+        double ratm2 = (is.romix - is.rorayl) * tgp2 + is.rorayl * tgp1;
+        double romeas1 = ratm1 + rsurf * tgtot;
+        double romeas2 = ratm2 + rsurf * tgtot;
+        double romeas3 = ratm3 + rsurf * tgtot;
 
 	/* computing integrated values over the spectral band */
         if (iwave.iwave == -2)
@@ -418,7 +418,7 @@ TransformInput compute()
         }
 
         
-	float alumeas = (float)(geom.xmus * swl * romeas2 / M_PI);
+	double alumeas = (double)(geom.xmus * swl * romeas2 / M_PI);
         fophsa = fophsa + is.phaa * coef;
         fophsr = fophsr + is.phar * coef;
         sasr = sasr + is.asray * coef;
@@ -474,15 +474,15 @@ TransformInput compute()
         seb = seb + coef;
 
 	/* 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);
-        float esn = tdir * dgtot;
-        float es = (float)(tdir * dgtot * geom.xmus * swl);
-        float ea0n = tdif * dgtot;
-        float ea0 = (float)(tdif * dgtot * geom.xmus * swl);
-        float ee0n = dgtot * avr * is.astot * is.dtott / (1 - avr * is.astot);
-        float ee0 = (float)(geom.xmus * swl * dgtot * avr * is.astot * is.dtott / (1 - avr * is.astot));
+        double tdir = (double)exp(-(is.tray + is.taer) / geom.xmus);
+        double tdif = is.dtott - tdir;
+        double etn = is.dtott * dgtot / (1 - avr * is.astot);
+        double esn = tdir * dgtot;
+        double es = (double)(tdir * dgtot * geom.xmus * swl);
+        double ea0n = tdif * dgtot;
+        double ea0 = (double)(tdif * dgtot * geom.xmus * swl);
+        double ee0n = dgtot * avr * is.astot * is.dtott / (1 - avr * is.astot);
+        double ee0 = (double)(geom.xmus * swl * dgtot * avr * is.astot * is.dtott / (1 - avr * is.astot));
 
         if (etn > accu3)
 	{
@@ -509,14 +509,14 @@ TransformInput compute()
 	}
 
 	/* output at satellite level */
-        float tmdir = (float)exp(-(is.tray + is.taer) / geom.xmuv);
-        float tmdif = is.utott - tmdir;
-        float xla0n = ratm2;
-        float xla0 = (float)(xla0n * geom.xmus * swl / M_PI);
-        float xltn = roc * is.dtott * tmdir * tgtot / (1 - avr * is.astot);
-        float xlt = (float)(xltn * geom.xmus * swl / M_PI);
-        float xlen = avr * is.dtott * tmdif * tgtot / (1 - avr * is.astot);
-        float xle = (float)(xlen * geom.xmus * swl / M_PI);
+        double tmdir = (double)exp(-(is.tray + is.taer) / geom.xmuv);
+        double tmdif = is.utott - tmdir;
+        double xla0n = ratm2;
+        double xla0 = (double)(xla0n * geom.xmus * swl / M_PI);
+        double xltn = roc * is.dtott * tmdir * tgtot / (1 - avr * is.astot);
+        double xlt = (double)(xltn * geom.xmus * swl / M_PI);
+        double xlen = avr * is.dtott * tmdif * tgtot / (1 - avr * is.astot);
+        double xle = (double)(xlen * geom.xmus * swl / M_PI);
         anr[0][0] = xla0n;
         anr[0][1] = xlen;
         anr[0][2] = xltn;
@@ -582,7 +582,7 @@ TransformInput compute()
     srotot = srotot / seb;
     alumet = alumet / sb;
     /*
-    float pizera = 0.0f;
+    double pizera = 0.0f;
     if(aero.iaer != 0) pizera = ssdaer / sodaer;
     */
     sodray = sodray / seb;

imagery/i.atcorr/6s.h

@@ -29,9 +29,9 @@ extern int init_6S(char* icnd_name);
   This requires lots of computations and therefore can be very
   time consuming.
 */
-extern void pre_compute_hv(const float height, const float vis);
-extern void pre_compute_v(const float vis);
-extern void pre_compute_h(const float height);
+extern void pre_compute_hv(const double height, const double vis);
+extern void pre_compute_v(const double vis);
+extern void pre_compute_h(const double height);
 
 struct TransformInput;
 /* Compute the input parameters used to do atmospheric correction on input values

imagery/i.atcorr/abstra.cpp

@@ -7,11 +7,11 @@
 #include "altitude.h"
 
 void
-wava6 (float a[8], const long int inu)
+wava6 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .011482f, .13183f,
+    static const double acr[2048] = { .011482f, .13183f,
 				     -.0038755f, 3.4491e-6f, -.0069899f, 9.3146e-6f, 15300.f, 15310.f,
 				     .0015124f, .19547f, .0028474f, -4.7616e-6f, .0017802f, -1.079e-5f,
 				     15310.f, 15320.f, .0092482f, .16207f, -.0025675f, 1.271e-5f, -.0027267f,
@@ -344,11 +344,11 @@ wava6 (float a[8], const long int inu)
 }	/* wava6 */
 
 void
-wava5 (float a[8], const long int inu)
+wava5 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 4.6416e-4f, .04653f,
+    static const double acr[2048] = { 4.6416e-4f, .04653f,
 				     .011484f, -5.0228e-5f, .0057564f, -2.8823e-5f, 12740.f, 12750.f,
 				     2.6026e-5f, .069686f, .0050381f, -3.0969e-5f, .0023565f, -2.6498e-5f,
 				     12750.f, 12760.f, 2.1016e-4f, .078469f, -.0024738f, -2.0423e-6f,
@@ -677,11 +677,11 @@ wava5 (float a[8], const long int inu)
 }	/* wava5 */
 
 void
-wava4 (float a[8], const long int inu)
+wava4 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .037011f, .34865f,
+    static const double acr[2048] = { .037011f, .34865f,
 				     .0071795f, -2.429e-5f, .0061217f, -2.5788e-5f, 10180.f, 10190.f,
 				     .096531f, .1963f, .0044353f, -2.7769e-5f, .0020496f, -1.902e-5f,
 				     10190.f, 10200.f, .11553f, .22356f, .0057418f, -2.861e-5f, .005252f,
@@ -1007,11 +1007,11 @@ wava4 (float a[8], const long int inu)
 }	/* wava4 */
 
 void
-wava3 (float a[8], const long int inu)
+wava3 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .092641f, .26739f,
+    static const double acr[2048] = { .092641f, .26739f,
 				     .0074828f, -3.6295e-5f, .0065918f, -3.6255e-5f, 7620.f, 7630.f, .24311f,
 				     .19859f, .0029686f, -1.983e-5f, .0023399f, -1.6807e-5f, 7630.f, 7640.f,
 				     .12025f, .11463f, .005982f, -3.2695e-5f, .00555f, -2.817e-5f, 7640.f,
@@ -1321,11 +1321,11 @@ wava3 (float a[8], const long int inu)
 }	/* wava3 */
 
 void
-wava2 (float a[8], const long int inu)
+wava2 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .32591f, .48473f,
+    static const double acr[2048] = { .32591f, .48473f,
 				     .010062f, 1.8245e-5f, .01189f, -1.2621e-5f, 5060.f, 5070.f, .73059f,
 				     .13181f, .010626f, 7.3795e-6f, .011376f, -1.7764e-5f, 5070.f, 5080.f,
 				     .39211f, .39522f, .01459f, -6.8376e-6f, .016326f, -3.165e-5f, 5080.f,
@@ -1637,11 +1637,11 @@ wava2 (float a[8], const long int inu)
 }	/* wava2 */
 
 void
-wava1 (float a[8], const long int inu)
+wava1 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 5.2155e-5f, .1088f,
+    static const double acr[2048] = { 5.2155e-5f, .1088f,
 				     .024708f, 5.6434e-5f, .028126f, -3.6504e-5f, 2500.f, 2510.f, 2.6024e-4f,
 				     .21216f, .025876f, 3.0026e-5f, .030504f, -6.2253e-5f, 2510.f, 2520.f,
 				     1.2221e-4f, .091374f, .023862f, -7.9891e-5f, .020651f, -8.5449e-5f,
@@ -1952,11 +1952,11 @@ wava1 (float a[8], const long int inu)
 }	/* wava1 */
 
 
-void dica3 (float a[8], const long int inu)
+void dica3 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 4.1135e-5f, .13491f,
+    static const double acr[2048] = { 4.1135e-5f, .13491f,
 				     .019511f, -8.8592e-5f, .017169f, -8.6383e-5f, 7620.f, 7630.f, 0.f, 0.f,
 				     0.f, 0.f, 0.f, 0.f, 7630.f, 7640.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7640.f,
 				     7650.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7650.f, 7660.f, 0.f, 0.f, 0.f, 0.f,
@@ -2256,11 +2256,11 @@ void dica3 (float a[8], const long int inu)
 }	/* dica3 */
 
 
-void dica2 (float a[8], const long int inu)
+void dica2 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .37011f, .18132f,
+    static const double acr[2048] = { .37011f, .18132f,
 				     .0098385f, -4.992e-5f, .0096965f, -3.9497e-5f, 5060.f, 5070.f, 1.7202f,
 				     .2316f, .0029954f, -2.2435e-5f, .0029757f, -9.2488e-6f, 5070.f, 5080.f,
 				     3.3606f, .25416f, -.0016977f, -4.0846e-6f, -.0013656f, 1.1658e-5f,
@@ -2562,12 +2562,12 @@ void dica2 (float a[8], const long int inu)
 }	/* dica2 */
 
 
-void dica1 (float a[8], const long int inu)
+void dica1 (double a[8], const long int inu)
 {
 
 
 
-    static const float acr[2048] = { 1.1446e-5f, .0020117f,
+    static const double acr[2048] = { 1.1446e-5f, .0020117f,
 				     -.0041334f, 3.2304e-6f, -.0069982f, 9.0084e-6f, 2500.f, 2510.f,
 				     1.9234e-5f, .0019311f, -.0017326f, -5.8646e-6f, -.0045311f,
 				     -6.0352e-7f, 2510.f, 2520.f, 9.202e-6f, .0017952f, .0034861f,
@@ -2873,11 +2873,11 @@ void dica1 (float a[8], const long int inu)
 
 
 void
-ozon1 (float a[8], const long int inu)
+ozon1 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .062007f, 2.4365f,
+    static const double acr[2048] = { .062007f, 2.4365f,
 				     -5.9503e-4f, -8.1198e-6f, -.0039418f, -2.4624e-6f, 2500.f, 2510.f,
 				     .023839f, 2.3534f, .0037377f, -6.15e-6f, .0015592f, -1.2727e-5f, 2510.f,
 				     2520.f, .0090127f, 1.2172f, -.0014733f, -4.7053e-6f, -.0042092f,
@@ -3172,11 +3172,11 @@ ozon1 (float a[8], const long int inu)
 
 
 void
-niox6 (float a[8], const long int inu)
+niox6 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 15300.f, 15310.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 15310.f, 15320.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 15320.f, 15330.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -3503,11 +3503,11 @@ niox6 (float a[8], const long int inu)
 }	/* niox6 */
 
 void
-niox5 (float a[8], const long int inu)
+niox5 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 12740.f, 12750.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 12750.f, 12760.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 12760.f, 12770.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -3834,11 +3834,11 @@ niox5 (float a[8], const long int inu)
 }	/* niox5 */
 
 void
-niox4 (float a[8], const long int inu)
+niox4 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 10180.f, 10190.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 10190.f, 10200.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 10200.f, 10210.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -4166,11 +4166,11 @@ niox4 (float a[8], const long int inu)
 }	/* niox4 */
 
 void
-niox3 (float a[8], const long int inu)
+niox3 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 7620.f, 7630.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7630.f, 7640.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 7640.f, 7650.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7650.f,
@@ -4461,11 +4461,11 @@ niox3 (float a[8], const long int inu)
 }	/* niox3 */
 
 void
-niox2 (float a[8], const long int inu)
+niox2 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { .072211f, .24584f,
+    static const double acr[2048] = { .072211f, .24584f,
 				     .0096738f, -5.1958e-5f, .0067533f, -4.7277e-5f, 5060.f, 5070.f, .21388f,
 				     .25456f, .0043318f, -3.1058e-5f, .0012217f, -2.5614e-5f, 5070.f, 5080.f,
 				     .57556f, .33263f, -2.6597e-4f, -1.2844e-5f, -.0033007f, -7.3238e-6f,
@@ -4751,11 +4751,11 @@ niox2 (float a[8], const long int inu)
 }	/* niox2 */
 
 void
-niox1 (float a[8], const long int inu)
+niox1 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 2.0198f, 1.2223f,
+    static const double acr[2048] = { 2.0198f, 1.2223f,
 				     .021725f, -7.4064e-5f, .021102f, -6.8716e-5f, 2500.f, 2510.f, 5.563f,
 				     .51358f, .018526f, -8.1387e-5f, .020173f, -7.5293e-5f, 2510.f, 2520.f,
 				     30.587f, .41845f, .010994f, -5.2858e-5f, .012658f, -4.4443e-5f, 2520.f,
@@ -5054,11 +5054,11 @@ niox1 (float a[8], const long int inu)
 
 
 void
-meth6 (float a[8], const long int inu)
+meth6 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 15300.f, 15310.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 15310.f, 15320.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 15320.f, 15330.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -5386,11 +5386,11 @@ meth6 (float a[8], const long int inu)
 }	/* meth6 */
 
 void
-meth5 (float a[8], const long int inu)
+meth5 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 12740.f, 12750.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 12750.f, 12760.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 12760.f, 12770.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -5715,11 +5715,11 @@ meth5 (float a[8], const long int inu)
 }	/* meth5 */
 
 void
-meth4 (float a[8], const long int inu)
+meth4 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 10180.f, 10190.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 10190.f, 10200.f,
 				     0.f,
 
@@ -6045,11 +6045,11 @@ meth4 (float a[8], const long int inu)
 }	/* meth4 */
 
 void
-meth3 (float a[8], const long int inu)
+meth3 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 7620.f, 7630.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7630.f, 7640.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 7640.f, 7650.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7650.f,
@@ -6339,11 +6339,11 @@ meth3 (float a[8], const long int inu)
 }	/* meth3 */
 
 void
-meth2 (float a[8], const long int inu)
+meth2 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 5060.f, 5070.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 5070.f, 5080.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 5080.f, 5090.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 5090.f,
@@ -6633,11 +6633,11 @@ meth2 (float a[8], const long int inu)
 }	/* meth2 */
 
 void
-meth1 (float a[8], const long int inu)
+meth1 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 1.4454f, .47807f,
+    static const double acr[2048] = { 1.4454f, .47807f,
 				     .0052823f, -3.0056e-5f, .002903f, -2.686e-5f, 2500.f, 2510.f, 8.7736f,
 				     .49348f, 3.8511e-4f, -6.0533e-6f, 1.0891e-4f, -9.3895e-6f, 2510.f,
 				     2520.f, 5.7188f, .51082f, 3.239e-4f, -7.2399e-6f, 1.6424e-4f,
@@ -6936,11 +6936,11 @@ meth1 (float a[8], const long int inu)
 
 
 void
-moca6 (float a[8], const long int inu)
+moca6 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 15300.f, 15310.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 15310.f, 15320.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 15320.f, 15330.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -7267,11 +7267,11 @@ moca6 (float a[8], const long int inu)
 }	/* moca6 */
 
 void
-moca5 (float a[8], const long int inu)
+moca5 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 12740.f, 12750.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 12750.f, 12760.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 12760.f, 12770.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -7597,9 +7597,9 @@ moca5 (float a[8], const long int inu)
 }	/* moca5 */
 
 
-void moca4 (float a[8], const long int inu)
+void moca4 (double a[8], const long int inu)
 {
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 10180.f, 10190.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 10190.f, 10200.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 10200.f, 10210.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -7927,11 +7927,11 @@ void moca4 (float a[8], const long int inu)
 }	/* moca4 */
 
 void
-moca3 (float a[8], const long int inu)
+moca3 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 7620.f, 7630.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7630.f, 7640.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 7640.f, 7650.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7650.f,
@@ -8233,11 +8233,11 @@ moca3 (float a[8], const long int inu)
 }	/* moca3 */
 
 void
-moca2 (float a[8], const long int inu)
+moca2 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 5060.f, 5070.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 5070.f, 5080.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 5080.f, 5090.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 5090.f,
@@ -8530,11 +8530,11 @@ moca2 (float a[8], const long int inu)
 }	/* moca2 */
 
 void
-moca1 (float a[8], const long int inu)
+moca1 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 2500.f, 2510.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 2510.f, 2520.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 2520.f, 2530.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 2530.f,
@@ -8840,11 +8840,11 @@ moca1 (float a[8], const long int inu)
 
 
 void
-oxyg6 (float a[8], const long int inu)
+oxyg6 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 15300.f, 15310.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 15310.f, 15320.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 15320.f, 15330.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -9167,9 +9167,9 @@ oxyg6 (float a[8], const long int inu)
     for(int i = 0; i < 8; i++) a[i] = acr[i + (inu << 3) - 8];
 }	/* oxyg6 */
 
-void oxyg5 (float a[8], const long int inu)
+void oxyg5 (double a[8], const long int inu)
 {
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 12740.f, 12750.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 12750.f, 12760.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 12760.f, 12770.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -9499,11 +9499,11 @@ void oxyg5 (float a[8], const long int inu)
 }	/* oxyg5 */
 
 void
-oxyg4 (float a[8], const long int inu)
+oxyg4 (double a[8], const long int inu)
 {
 
 
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 10180.f, 10190.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 10190.f, 10200.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 0.f, 10200.f, 10210.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
@@ -9828,9 +9828,9 @@ oxyg4 (float a[8], const long int inu)
 }	/* oxyg4 */
 
 void
-oxyg3 (float a[8], const long int inu)
+oxyg3 (double a[8], const long int inu)
 {
-    static const float acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
+    static const double acr[2048] = { 0.f, 0.f, 0.f, 0.f, 0.f,
 				     0.f, 7620.f, 7630.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7630.f, 7640.f, 0.f,
 				     0.f,
 				     0.f, 0.f, 0.f, 0.f, 7640.f, 7650.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 7650.f,
@@ -10146,28 +10146,28 @@ oxyg3 (float a[8], const long int inu)
 
 
 void abstra (const AtmosModel& atms, const Altitude& alt,
-	     const float wl, const float xmus, const float xmuv,
-	     const float uw, const float uo3, float& uwus, float& uo3us,
-	     const float uwpl, const float uo3pl, const float uwusp,
-	     const float uo3usp, AbstraStruct& as )
+	     const double wl, const double xmus, const double xmuv,
+	     const double uw, const double uo3, double& uwus, double& uo3us,
+	     const double uwpl, const double uo3pl, const double uwusp,
+	     const double uo3usp, AbstraStruct& as )
 { 
     // transmittance calculation for ozone, water vapor,
     // carbon dioxyde and oxygen.
 	 
-    float tnu[10][3],a[8],rm[34],r2[34],r3[34],tp[34],rat[10];
-    float rmpl[34],r2pl[34],r3pl[34],ratpl[10];
-    float dtcont,utcont,ttcont;
-    float v,te2,phi,psi,uu,u,up,uud,uut,uuu;
-    float ud,ut,upd,upt,udp,updp,udtp,updtp;
-    float ds2,uupl,upl,uppl,ah2o;
-    float xh,xi,xd,ako3,test1,test2,test3,udt,atest;
-    float updt,tt,y,utt,uptt,tn;
-    float ds,te,roair;
+    double tnu[10][3],a[8],rm[34],r2[34],r3[34],tp[34],rat[10];
+    double rmpl[34],r2pl[34],r3pl[34],ratpl[10];
+    double dtcont,utcont,ttcont;
+    double v,te2,phi,psi,uu,u,up,uud,uut,uuu;
+    double ud,ut,upd,upt,udp,updp,udtp,updtp;
+    double ds2,uupl,upl,uppl,ah2o;
+    double xh,xi,xd,ako3,test1,test2,test3,udt,atest;
+    double updt,tt,y,utt,uptt,tn;
+    double ds,te,roair;
     double ptest, ptest1;
     int iv,id,idgaz,inu = 0,n,nh;
 
     int ivli[6] = { 2500,5060,7620,10180,12740,15300 };
-    float co3[102] = {
+    double co3[102] = {
 	4.50e-03, 8.00e-03, 1.07e-02, 1.10e-02, 1.27e-02, 1.71e-02,
 	2.00e-02, 2.45e-02, 3.07e-02, 3.84e-02, 4.78e-02, 5.67e-02,
 	6.54e-02, 7.62e-02, 9.15e-02, 1.00e-01, 1.09e-01, 1.20e-01,
@@ -10187,7 +10187,7 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	1.57e+01, 1.20e+01, 1.00e+01, 8.80e+00, 8.30e+00, 8.60e+00
     };
 
-    float cch2o[15] = {
+    double cch2o[15] = {
 	0.00,0.19,0.15,0.12,0.10,
 	0.09,0.10,0.12,0.15,0.17,
 	0.20,0.24,0.28,0.33,0.00
@@ -10236,22 +10236,22 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
     }
 		
     /* constants determination */
-    const float p0 = 1013.25f;
-    const float g = 98.1f;
+    const double p0 = 1013.25f;
+    const double g = 98.1f;
 
-    const float t0 = 250.f;
+    const double t0 = 250.f;
  
     /* volumic mass in kilogrammes per m3 */
     ds = 0;
     te = 0;
     roair = 0;
-    const float air = 0.028964 / 0.0224;
-    const float roco2 = 0.044 / 0.0224;
-    const float rmo2 = 0.032 / 0.0224;
-    const float rmo3 = 0.048 / 0.0224;
-    const float rmn2o = 0.044 / 0.0224;
-    const float rmch4 = 0.016 / 0.0224;
-    const float rmco  = 0.028 / 0.0224;
+    const double air = 0.028964 / 0.0224;
+    const double roco2 = 0.044 / 0.0224;
+    const double rmo2 = 0.032 / 0.0224;
+    const double rmo3 = 0.048 / 0.0224;
+    const double rmn2o = 0.044 / 0.0224;
+    const double rmch4 = 0.016 / 0.0224;
+    const double rmco  = 0.028 / 0.0224;
 	
     uwus = 1.424;
     uo3us = .344;
@@ -10270,7 +10270,7 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	rat[9] = uw / uwus;
     }
  
-    v = (float)(1e+04 / wl);
+    v = (double)(1e+04 / wl);
     iv = (int)(v / 5);
     iv = iv * 5;
     id = ((iv - 2500) / 10) / 256 + 1;
@@ -10361,8 +10361,8 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	    tp[k] = (atms.t[k] + atms.t[k + 1]) / 2.f;
 	    te = tp[k] - t0;
 	    te2 = te * te;
-	    phi = (float)exp(a[2] * te + a[3] * te2);
-	    psi = (float)exp(a[4] * te + a[5] * te2);
+	    phi = (double)exp(a[2] * te + a[3] * te2);
+	    psi = (double)exp(a[4] * te + a[5] * te2);
 	    if(idgaz == 1) rm[k] = atms.wh[k] / (roair * 1000);
 	    if(idgaz == 2) rm[k] = 3.3e-04f * roco2 / air;
 	    if(idgaz == 3) rm[k] = 0.20947f * rmo2 / air;
@@ -10424,8 +10424,8 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 		tp[k]  =  (alt.plane_sim.tpl[k] + alt.plane_sim.tpl[k + 1]) / 2;
 		te = tp[k] - t0;
 		te2 = te * te;
-		phi = (float)exp(a[2] * te + a[3] * te2);
-		psi = (float)exp(a[4] * te + a[5] * te2);
+		phi = (double)exp(a[2] * te + a[3] * te2);
+		psi = (double)exp(a[4] * te + a[5] * te2);
 		if(idgaz == 1) rmpl[k] = alt.plane_sim.whpl[k] / (roair * 1000);
 		if(idgaz == 2) rmpl[k] = 3.3e-04f * roco2 / air;
 		if(idgaz == 3) rmpl[k] = 0.20947f * rmo2 / air;
@@ -10509,11 +10509,11 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	    {
 		xi = (v - 2350) / 50 + 1;
 		nh = (int)(xi + 1.001f);
-		xh = xi - float(nh);
+		xh = xi - double(nh);
 		ah2o = cch2o[nh-1] + xh * (cch2o[nh-1]-cch2o[nh-2]);
-		dtcont = (float)exp(-ah2o * uud);
-		utcont = (float)exp(-ah2o * uuu);
-		ttcont = (float)exp(-ah2o * uut);
+		dtcont = (double)exp(-ah2o * uud);
+		utcont = (double)exp(-ah2o * uuu);
+		ttcont = (double)exp(-ah2o * uut);
 	    }
 
 	    if (!((idgaz == 1) || (iv < 13000))) 
@@ -10529,7 +10529,7 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 		}
 
 		n = (int)(xi + 1.001);
-		xd = xi-float(n);
+		xd = xi-double(n);
 		ako3 = co3[n-1] + xd * (co3[n-1] - co3[n-2]);
 		test1 = ako3 * uud;
 		test2 = ako3 * uuu;
@@ -10541,9 +10541,9 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 		if(test2 > 86.0) test2 = 86.0;
 		if(test3 > 86.0) test3 = 86.0;
 	 
-		tnu[3][0] = (float)exp(-test1);
-		tnu[3][1] = (float)exp(-test2);
-		tnu[3][2] = (float)exp(-test3);
+		tnu[3][0] = (double)exp(-test1);
+		tnu[3][1] = (double)exp(-test2);
+		tnu[3][2] = (double)exp(-test3);
 
 		continue;
 	    }
@@ -10567,9 +10567,9 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	updt = upd;
 	if(ud == 0 && upd == 0.) updt = 1;
 	tt = 1 + 4 * (a[0] / atest) * ((ud * ud) / updt);
-	y = (float)(-tn * (sqrt(tt) - 1));
-	if(idgaz == 1) y = (float)(-a[0] * ud / sqrt(1 + (a[0] / atest) * (ud * ud / updt)));
-	tnu[idgaz-1][0] = (float)exp(y);
+	y = (double)(-tn * (sqrt(tt) - 1));
+	if(idgaz == 1) y = (double)(-a[0] * ud / sqrt(1 + (a[0] / atest) * (ud * ud / updt)));
+	tnu[idgaz-1][0] = (double)exp(y);
 		
 			
 	/* upward path modified to take account for plane content */
@@ -10583,9 +10583,9 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	updtp = updp;
 	if(udp == 0 && updp == 0.) updtp = 1;
 	tt = 1 + 4 * (a[0] / atest) * ((udp * udp) / updtp);
-	y = (float)(-tn * (sqrt(tt) - 1));
-	if(idgaz == 1) y = (float)(-a[0] * udp / sqrt(1 + (a[0] / atest) * (udp * udp / updtp)));
-	tnu[idgaz-1][1] = (float)exp(y);
+	y = (double)(-tn * (sqrt(tt) - 1));
+	if(idgaz == 1) y = (double)(-a[0] * udp / sqrt(1 + (a[0] / atest) * (udp * udp / updtp)));
+	tnu[idgaz-1][1] = (double)exp(y);
 
 	/* total(down + up) path modified on the way up */
 	ut = u / xmus + upl / xmuv;
@@ -10596,26 +10596,26 @@ void abstra (const AtmosModel& atms, const Altitude& alt,
 	uptt = upt;
 	if(ut == 0 && upt == 0.) uptt = 1;
 	tt = 1 + 4 * (a[0] / atest) * ((ut * ut) / uptt);
-	y = (float)(-tn * (sqrt(tt) - 1));
-	if(idgaz == 1) y = (float)(-a[0] * ut / sqrt(1 + (a[0] / atest) * (ut * ut / uptt)));
-	tnu[idgaz-1][2] = (float)exp(y);
+	y = (double)(-tn * (sqrt(tt) - 1));
+	if(idgaz == 1) y = (double)(-a[0] * ut / sqrt(1 + (a[0] / atest) * (ut * ut / uptt)));
+	tnu[idgaz-1][2] = (double)exp(y);
     }   
 
     ptest1 = tnu[0][0] * dtcont;
     ptest = ptest1;
-    if (ptest > 1e-10) as.dtwava = (float)ptest;
+    if (ptest > 1e-10) as.dtwava = (double)ptest;
     else as.dtwava = 0;
 	
     ptest1 = tnu[0][1] * utcont;
     ptest = ptest1;
 
-    if (ptest > 1e-10) as.utwava = (float)ptest;
+    if (ptest > 1e-10) as.utwava = (double)ptest;
     else as.utwava = 0;
 	
     ptest1 = tnu[0][2] * ttcont;
     ptest = ptest1;
 
-    if (ptest > 1e-10) as.ttwava = (float)ptest;
+    if (ptest > 1e-10) as.ttwava = (double)ptest;
     else as.ttwava = 0;
 	
     as.dtdica = tnu[1][0];

imagery/i.atcorr/abstra.h

@@ -3,27 +3,27 @@
 
 struct AbstraStruct
 {
-	float dtwava; /* downward absorption water vapor dtwava */
-	float dtozon; /* downward absorption ozone       dtozon */
-	float dtdica; /* downward absorption carbon diox dtdica */
-	float dtoxyg; /* downward absorption oxygen      dtoxyg */
-	float dtniox; /* downward absorption nitrous oxi dtniox */
-	float dtmeth; /* downward absorption methane     dtmeth */
-	float dtmoca; /* downward absorption carbon mono dtmoca */
-	float utwava; /* upward absorption water vapor   utwava */
-	float utozon; /* upward absorption ozone         utozon */
-	float utdica; /* upward absorption carbon diox   utdica */
-	float utoxyg; /* upward absorption oxygen        utoxyg */
-	float utniox; /* upward   absorption nitrous oxi utniox */
-	float utmeth; /* upward   absorption methane     utmeth */
-	float utmoca; /* upward   absorption carbon mono utmoca */
-	float ttwava; /* total(on the two paths ) absorption water vapor ttwava */
-	float ttozon; /* total(on the two paths ) absorption ozone       ttozon */
-	float ttdica; /* total(on the two paths ) absorption carbon diox ttdica */
-	float ttoxyg; /* total(on the two paths ) absorption oxygen      ttoxyg */
-	float ttniox; /* total    absorption nitrous oxi ttniox */
-	float ttmeth; /* total    absorption methane     ttmeth */
-	float ttmoca; /* total    absorption carbon mono ttmoca */
+	double dtwava; /* downward absorption water vapor dtwava */
+	double dtozon; /* downward absorption ozone       dtozon */
+	double dtdica; /* downward absorption carbon diox dtdica */
+	double dtoxyg; /* downward absorption oxygen      dtoxyg */
+	double dtniox; /* downward absorption nitrous oxi dtniox */
+	double dtmeth; /* downward absorption methane     dtmeth */
+	double dtmoca; /* downward absorption carbon mono dtmoca */
+	double utwava; /* upward absorption water vapor   utwava */
+	double utozon; /* upward absorption ozone         utozon */
+	double utdica; /* upward absorption carbon diox   utdica */
+	double utoxyg; /* upward absorption oxygen        utoxyg */
+	double utniox; /* upward   absorption nitrous oxi utniox */
+	double utmeth; /* upward   absorption methane     utmeth */
+	double utmoca; /* upward   absorption carbon mono utmoca */
+	double ttwava; /* total(on the two paths ) absorption water vapor ttwava */
+	double ttozon; /* total(on the two paths ) absorption ozone       ttozon */
+	double ttdica; /* total(on the two paths ) absorption carbon diox ttdica */
+	double ttoxyg; /* total(on the two paths ) absorption oxygen      ttoxyg */
+	double ttniox; /* total    absorption nitrous oxi ttniox */
+	double ttmeth; /* total    absorption methane     ttmeth */
+	double ttmoca; /* total    absorption carbon mono ttmoca */
 };
 
 struct AtmosModel;
@@ -36,9 +36,9 @@ The total transmission is put equal to the simple product of each ones. The spec
 equal to 10 cm-1.
 iinf*/
 void abstra (const AtmosModel& atms, const Altitude& alt,
-			 const float wl, const float xmus, const float xmuv,
-			 const float uw, const float uo3, float& uwus, float& uo3us,
-			 const float uwpl, const float uo3pl, const float uwusp,
-			 const float uo3usp, AbstraStruct& as );
+			 const double wl, const double xmus, const double xmuv,
+			 const double uw, const double uo3, double& uwus, double& uo3us,
+			 const double uwpl, const double uo3pl, const double uwusp,
+			 const double uo3usp, AbstraStruct& as );
 
 #endif /* ABSTRA_H */

imagery/i.atcorr/aerosolconcentration.cpp

@@ -32,17 +32,17 @@ void AerosolConcentration::parse(const long int _iaer, const AtmosModel& atms)
     {
 	cin >> taer55;
 	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
-	v = (float)(exp(-log(taer55/2.7628f)/0.79902f));
+	v = (double)(exp(-log(taer55/2.7628f)/0.79902f));
     }
     else if(v > 0) oda550(v, atms);
 }
 
-void AerosolConcentration::oda550(const float vis, const AtmosModel& atms)
+void AerosolConcentration::oda550(const double vis, 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] = {
+    static const double 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,
@@ -55,7 +55,7 @@ void AerosolConcentration::oda550(const float vis, const AtmosModel& atms)
 
     /* vertical repartition of aerosol density for v=5km */
     /*                ( in nbr of part/cm3 ) */
-    static const float an5[34] = {
+    static const double 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,
@@ -71,17 +71,17 @@ void AerosolConcentration::oda550(const float vis, const AtmosModel& atms)
 
     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]);
+	double dz = atms.z[k+1] - atms.z[k];
+	double az = (115.f / 18.f) * (an5[k] - an23[k]);
+	double 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);
+	double bz = (5.f * an5[k] / 18.f) - (23.f * an23[k] / 18.f);
+	double bz1 = (5.f * an5[k+1] / 18.f) - (23.f * an23[k+1] / 18.f);
 
-	float bnz = az / vis - bz;
-	float bnz1 = az1 / vis - bz1;
+	double bnz = az / vis - bz;
+	double bnz1 = az1 / vis - bz1;
 
-	float ev = (float)(dz * exp((log(bnz) + log(bnz1)) / 2));
+	double ev = (double)(dz * exp((log(bnz) + log(bnz1)) / 2));
 	taer55 += ev * sigma * 1.0e-03f;
     }
 }

imagery/i.atcorr/aerosolconcentration.h

@@ -32,17 +32,17 @@ struct AtmosModel;
 struct AerosolConcentration
 {
 	/* aerosol concentration parameters */
-    float taer55;
+    double taer55;
 
 private:
     long int iaer;
-    float v;
+    double v;
     void parse(const long int iaer, const AtmosModel &atms);
-    void oda550(const float v, const AtmosModel &atms);
+    void oda550(const double v, const AtmosModel &atms);
 
 public:
     /* Set the visibility, this will overide any previous estimates of taer55 */
-    void set_visibility (const float vis, const AtmosModel &atms) { if(vis > 0) oda550(vis, atms); }
+    void set_visibility (const double vis, const AtmosModel &atms) { if(vis > 0) oda550(vis, atms); }
     void print();
     static AerosolConcentration Parse(const long int iaer, const AtmosModel &atms);
 };

imagery/i.atcorr/aerosolmodel.cpp

@@ -7,7 +7,7 @@ extern "C" {
 #include "aerosolmodel.h"
 #include "atmosmodel.h"
 #ifdef WIN32
-#pragma warning(disable:4305)	/* disable warning about initialization of a float by a double */
+#pragma warning(disable:4305)	/* disable warning about initialization of a double by a double */
 #endif /* WIN32 */
 
 /* (background desert model...) */
@@ -57,7 +57,7 @@ void AerosolModel::soot()
   isotropic sphere, the physical properties of particles whose sizes are 
   comparable to or larger than the wavelength, and to generate mixture of 
   dry particles. */
-void AerosolModel::mie(float (&ex)[4][10], float (&sc)[4][10], float (&asy)[4][10])
+void AerosolModel::mie(double (&ex)[4][10], double (&sc)[4][10], double (&asy)[4][10])
 {
     double np[4];
     double ext[10][4];
@@ -180,8 +180,8 @@ void AerosolModel::mie(float (&ex)[4][10], float (&sc)[4][10], float (&asy)[4][1
 	{
 	    ext[j][i] /= np[i] * 1000;
 	    sca2[j][i] /= np[i] * 1000;
-	    ex[0][j] += (float)(mie_in.cij[i] * ext[j][i]);
-	    sc[0][j] += (float)(mie_in.cij[i] * sca2[j][i]);
+	    ex[0][j] += (double)(mie_in.cij[i] * ext[j][i]);
+	    sc[0][j] += (double)(mie_in.cij[i] * sca2[j][i]);
 	}
 
     /* computation of the phase function and the asymetry coefficient
@@ -196,14 +196,14 @@ void AerosolModel::mie(float (&ex)[4][10], float (&sc)[4][10], float (&asy)[4][1
 	{
 	    sixs_aerbas.usr_ph[j][k] = 0;
 	    for(i = 0; i < mie_in.icp; i++)
-		sixs_aerbas.usr_ph[j][k] += (float)(mie_in.cij[i] * p1[j][i][k] / np[i] / 1000);
+		sixs_aerbas.usr_ph[j][k] += (double)(mie_in.cij[i] * p1[j][i][k] / np[i] / 1000);
 		
-	    sixs_aerbas.usr_ph[j][k] += (float)sc[0][j];
+	    sixs_aerbas.usr_ph[j][k] += (double)sc[0][j];
 	    asy_n += sixs_sos.cgaus[k] * sixs_aerbas.usr_ph[j][k] * sixs_sos.pdgs[k] / 10.;
 	    asy_d += sixs_aerbas.usr_ph[j][k] * sixs_sos.pdgs[k] / 10.;
 	}
 
-	asy[0][j] = (float)(asy_n / asy_d);
+	asy[0][j] = (double)(asy_n / asy_d);
     }
 
     sixs_aerbas.ph = &sixs_aerbas.usr_ph;
@@ -499,7 +499,7 @@ void AerosolModel::save()
   by the user (see SUBROUTINES MIE and EXSCPHASE).
   These models don't correspond to a mixture of the four basic components.
 */
-void AerosolModel::aeroso(const float xmud)
+void AerosolModel::aeroso(const double xmud)
 {
 /* sra basic components for aerosol model, extinction coefficients are */
 /* in km-1. */
@@ -511,7 +511,7 @@ void AerosolModel::aeroso(const float xmud)
     static const double ni[4] = { 54.734, 1868550., 276.05, 1805820. };
 
     /* i: 1=dust-like 2=water-soluble 3=oceanic 4=soot */
-    static const float s_ex[4][10] =
+    static const double s_ex[4][10] =
 	{
 	    {0.1796674e-01,0.1815135e-01,0.1820247e-01,0.1827016e-01,0.1842182e-01,
 	     0.1853081e-01,0.1881427e-01,0.1974608e-01,0.1910712e-01,0.1876025e-01},
@@ -523,7 +523,7 @@ void AerosolModel::aeroso(const float xmud)
 	     0.3966041e-06,0.2965532e-06,0.1493927e-06,0.1017134e-06,0.6065031e-07}
 	};
 
-    static const float s_sc[4][10] =
+    static const double s_sc[4][10] =
 	{
 	    {0.1126647e-01,0.1168918e-01,0.1180978e-01,0.1196792e-01,0.1232056e-01,
 	     0.1256952e-01,0.1319347e-01,0.1520712e-01,0.1531952e-01,0.1546761e-01},
@@ -535,44 +535,44 @@ void AerosolModel::aeroso(const float xmud)
 	     0.6469735e-07,0.3610638e-07,0.6227224e-08,0.1779378e-08,0.3050002e-09}
 	};
  
-    static const float ex2[10] = 
+    static const double ex2[10] = 
 	{ 
 	    43.83631f, 42.12415f, 41.57425f, 40.85399f, 39.1404f, 
 	    37.89763f, 34.67506f, 24.59f, 17.96726f, 10.57569f
 	};
 
-    static const float sc2[10] = 
+    static const double sc2[10] = 
 	{ 
 	    40.28625f, 39.04473f, 38.6147f, 38.03645f, 36.61054f, 
 	    35.54456f, 32.69951f, 23.41019f, 17.15375f,10.09731f
 	};
 
-    static const float ex3[10] = 
+    static const double ex3[10] = 
 	{ 
 	    95397.86f, 75303.6f, 70210.64f, 64218.28f, 52430.56f, 
 	    45577.68f, 31937.77f, 9637.68f, 3610.691f, 810.5614f
 	};
 
-    static const float sc3[10] = 
+    static const double sc3[10] = 
 	{ 
 	    92977.9f, 73397.17f, 68425.49f,	62571.8f, 51049.87f, 
 	    44348.77f, 31006.21f, 9202.678f, 3344.476f,	664.1915f
 	};
   
-    static const float ex4[10] = 
+    static const double ex4[10] = 
 	{ 
 	    54273040.f, 61981440.f, 63024320.f, 63489470.f, 61467600.f,
 	    58179720.f, 46689090.f, 15190620.f, 5133055.f, 899859.4f
 	};
 
   
-    static const float sc4[10] = 
+    static const double sc4[10] = 
 	{ 
 	    54273040.f, 61981440.f, 63024320.f, 63489470.f, 61467600.f, 
 	    58179720.f, 46689090.f, 15190620.f, 5133055.f, 899859.4f
 	};
 	
-    static const float s_asy[4][10] = 
+    static const double s_asy[4][10] = 
 	{
 	    {0.896,0.885,0.880,0.877,0.867,0.860,0.845,0.836,0.905,0.871},
 	    {0.642,0.633,0.631,0.628,0.621,0.616,0.610,0.572,0.562,0.495},
@@ -580,21 +580,21 @@ void AerosolModel::aeroso(const float xmud)
 	    {0.397,0.359,0.348,0.337,0.311,0.294,0.253,0.154,0.103,0.055}
 	};
 
-    static const float asy2[10] = { .718f, .712f, .71f, .708f, .704f, .702f, .696f, .68f, .668f, .649f };
+    static const double asy2[10] = { .718f, .712f, .71f, .708f, .704f, .702f, .696f, .68f, .668f, .649f };
 
-    static const float asy3[10] = { .704f, .69f, .686f, .68f, .667f, .659f, .637f, .541f, .437f, .241f };
-    static const float asy4[10] = { .705f, .744f, .751f, .757f, .762f, .759f, .737f, .586f, .372f, .139f };
+    static const double asy3[10] = { .704f, .69f, .686f, .68f, .667f, .659f, .637f, .541f, .437f, .241f };
+    static const double asy4[10] = { .705f, .744f, .751f, .757f, .762f, .759f, .737f, .586f, .372f, .139f };
 
     /* local */
     double coef;
-    float sigm;
+    double sigm;
     double sumni;
     double dd[4][10];
     double pha[5][10][83];
 
-    float ex[4][10];
-    float sc[4][10];
-    float asy[4][10];
+    double ex[4][10];
+    double sc[4][10];
+    double asy[4][10];
 
     int i;	/* crappy VS6 */
     /* initialize ex, sc & asy */
@@ -636,7 +636,7 @@ void AerosolModel::aeroso(const float xmud)
     {
 	load();
 	for(i = 0; i < 10; i++) 
-	    sixs_aer.phase[i] = (float)(sixs_sos.phasel[i][j1] + 
+	    sixs_aer.phase[i] = (double)(sixs_sos.phasel[i][j1] + 
 					coef*(sixs_sos.phasel[i][j1]-sixs_sos.phasel[i][j1+1]));
 	return;
     }
@@ -745,11 +745,11 @@ void AerosolModel::aeroso(const float xmud)
 	sumni = 0.f;
 	sigm = 0.f;
 
-	for(i = 0; i < 4; i++) sigm+=(float)(c[i]/vi[i]);
+	for(i = 0; i < 4; i++) sigm+=(double)(c[i]/vi[i]);
 
 	/* cij coefficients calculation */
 	for(i = 0; i < 4; i++) {
-	    mie_in.cij[i] = (float)(c[i]/vi[i]/sigm);
+	    mie_in.cij[i] = (double)(c[i]/vi[i]/sigm);
 	    sumni += mie_in.cij[i]/ni[i];
 	}
 
@@ -762,13 +762,13 @@ void AerosolModel::aeroso(const float xmud)
     {
 	for(int j = 0; j < mie_in.icp; j++)
 	{
-	    sixs_aer.ext[i] +=		(float)(ex[j][i] * mie_in.cij[j]);
-	    sca[i] +=				(float)(sc[j][i] * mie_in.cij[j]);
-	    sixs_aer.gasym[i] +=	(float)(sc[j][i] * mie_in.cij[j] * asy[j][i]);
-	    sixs_aer.phase[i] +=	(float)(sc[j][i] * mie_in.cij[j] * dd[j][i]);
+	    sixs_aer.ext[i] +=		(double)(ex[j][i] * mie_in.cij[j]);
+	    sca[i] +=				(double)(sc[j][i] * mie_in.cij[j]);
+	    sixs_aer.gasym[i] +=	(double)(sc[j][i] * mie_in.cij[j] * asy[j][i]);
+	    sixs_aer.phase[i] +=	(double)(sc[j][i] * mie_in.cij[j] * dd[j][i]);
 
 	    for(int k = 0; k < 83; k++)
-		sixs_sos.phasel[i][k] += (float)(sc[j][i] * mie_in.cij[j] * pha[j][i][k]);
+		sixs_sos.phasel[i][k] += (double)(sc[j][i] * mie_in.cij[j] * pha[j][i][k]);
 	}
 
 	sixs_aer.ome[i] = sca[i]/sixs_aer.ext[i];
@@ -777,14 +777,14 @@ void AerosolModel::aeroso(const float xmud)
 
 	for(int k = 0; k < 83; k++)	sixs_sos.phasel[i][k] /= sca[i];
 
-	sixs_aer.ext[i] *= (float)nis;
-	sca[i] *= (float)nis;
+	sixs_aer.ext[i] *= (double)nis;
+	sca[i] *= (double)nis;
     }
 
     if (filename.size() != 0 && iaer >= 8 && iaer <= 11) save();
 }
 
-void AerosolModel::parse(const float xmud)
+void AerosolModel::parse(const double xmud)
 {
     cin >> iaer;
     cin.ignore(numeric_limits<int>::max(),'\n');
@@ -937,7 +937,7 @@ void AerosolModel::parse(const float xmud)
 
 	    double sq = mie_in.rsunph[i]*mie_in.rsunph[i];
 	    const double ln10 = 2.3025850929940456840179914546844;
-	    mie_in.nrsunph[i] = (float)(mie_in.nrsunph[i]/(sq*sq)/ln10);
+	    mie_in.nrsunph[i] = (double)(mie_in.nrsunph[i]/(sq*sq)/ln10);
 	}
 	mie_in.rmin=mie_in.rsunph[0];
 	mie_in.rmax=mie_in.rsunph[mie_in.irsunph-1]+1e-07f;
@@ -1159,7 +1159,7 @@ void AerosolModel::print()
     }
 }
 
-AerosolModel AerosolModel::Parse(const float xmud)
+AerosolModel AerosolModel::Parse(const double xmud)
 {
     AerosolModel aero;
     aero.parse(xmud);

imagery/i.atcorr/aerosolmodel.h

@@ -85,15 +85,15 @@
 struct AerosolModel
 {
 	long int iaer;	/* aerosol model */
-	float c[4];
+	double c[4];
 
 private:
 	double nis;
-	float sca[10];
+	double sca[10];
 	long int iaerp;
 
 	/* methods */
-	void aeroso(const float xmud);
+	void aeroso(const double xmud);
 
 	string filename;
 	void load();
@@ -110,34 +110,34 @@ private:
 
 	struct Mie_in
 	{
-		float rmax;
-		float rmin;
-		float rn[10][4];
-		float ri[10][4];
-		float x1[4];
-		float x2[4];
-		float x3[4];
-		float cij[4];
-		float rsunph[50];
-		float nrsunph[50];
+		double rmax;
+		double rmin;
+		double rn[10][4];
+		double ri[10][4];
+		double x1[4];
+		double x2[4];
+		double x3[4];
+		double cij[4];
+		double rsunph[50];
+		double nrsunph[50];
 
 		long int icp;
 		long int irsunph;
 	};
 
 	Mie_in mie_in;
-	void mie(float (&ex)[4][10], float (&sc)[4][10], float (&asy)[4][10]);
+	void mie(double (&ex)[4][10], double (&sc)[4][10], double (&asy)[4][10]);
 	void exscphase(const double alpha, const double nr, 
 				   const double ni, double& Qext, 
 				   double& Qsca, double (&p11)[83]);
 
-	void parse(const float xmud);
+	void parse(const double xmud);
 
 	/* format 132 */
 	void print132(string s);
 public:
 	void print();
-	static AerosolModel Parse(const float xmud);
+	static AerosolModel Parse(const double xmud);
 };
 
 

imagery/i.atcorr/altitude.cpp

@@ -14,7 +14,7 @@
    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)
+void Altitude::pressure(AtmosModel& atms, double& uw, double& uo3)
 {
     /* log linear interpolation */
     if(xps >= 100) xps = 99.99f;
@@ -25,17 +25,17 @@ void Altitude::pressure(AtmosModel& atms, float& uw, float& uo3)
     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);
+    double xa = (double)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
+    double xb = (double)(atms.z[isup] - xa * log(atms.p[isup]));
+    double ps = (double)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]);
+    double xalt = xps;
+    double 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]);
+    double 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]);
+    double 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
@@ -68,14 +68,14 @@ void Altitude::pressure(AtmosModel& atms, float& uw, float& uo3)
     /* 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;
+    const double g = 98.1f;
+    const double air = 0.028964f/0.0224f;
+    const double ro3 = 0.048f/0.0224f;
 
-    float rmwh[34];
-    float rmo3[34];
+    double rmwh[34];
+    double rmo3[34];
     int k;
-    float roair, ds;
+    double roair, ds;
 
     k = 0;
     roair = air * 273.16f * atms.p[k] / (atms.t[k] * 1013.25f);
@@ -122,17 +122,17 @@ void Altitude::presplane(AtmosModel& atms)
     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));
+    double xa = (double)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
+    double xb = (double)(atms.z[isup] - xa * log(atms.p[isup]));
+    double ps = (double)(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]);
+    double xalt = xpp;
+    double 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]);
+    double 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]);
+    double 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
@@ -161,18 +161,18 @@ void Altitude::presplane(AtmosModel& atms)
        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];
+    const double g = 98.1f;
+    const double air = 0.028964f/0.0224f;
+    const double ro3 = 0.048f/0.0224f;
+    double rt = 0;
+    double rp = 0;
+
+    double rmo3[34];
+    double 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]));
+	double roair = (double)(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]);
@@ -183,7 +183,7 @@ void Altitude::presplane(AtmosModel& atms)
     ftray = rp / rt;
     for(k = 1; k < 33; k++)
     {
-	float ds = (plane_sim.ppl[k-1] - plane_sim.ppl[k]) / plane_sim.ppl[0];
+	double 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;
     }
@@ -198,8 +198,8 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
     xps = original_xps;
     xpp = original_xpp;
 
-    float uwus;
-    float uo3us;
+    double uwus;
+    double uo3us;
     if(xps <= 0)
     {
 	xps = 0;
@@ -265,7 +265,7 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
 	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)));
+	    taer55p = (double)(aerocon.taer55 * (1 - exp(-palt / 2)));
 	}
 	else
 	{
@@ -273,10 +273,10 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
 	    double sham = exp(-palt / 4);
 	    double sha = 1 - (taer55p / aerocon.taer55);
 
-	    if( sha >= sham) taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
+	    if( sha >= sham) taer55p = (double)(aerocon.taer55 * (1 - exp(-palt / 4)));
 	    else {
 		sha = -palt/log(sha);
-		taer55p = (float)(aerocon.taer55 * (1 - exp(-palt/sha)));
+		taer55p = (double)(aerocon.taer55 * (1 - exp(-palt/sha)));
 	    }
 	}
     }
@@ -287,8 +287,8 @@ void Altitude::update_hv(AtmosModel & atms, const AerosolConcentration & aerocon
     xps = original_xps;
     xpp = original_xpp;
 
-    float uwus;
-    float uo3us;
+    double uwus;
+    double uo3us;
 
     if (xps <= 0) {
 	xps = 0;
@@ -346,7 +346,7 @@ void Altitude::update_hv(AtmosModel & atms, const AerosolConcentration & aerocon
 
 	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)));
+	    taer55p = (double)(aerocon.taer55 * (1 - exp(-palt / 2)));
 	}
 	else {
 	    /* compute effective scale heigh */
@@ -354,10 +354,10 @@ void Altitude::update_hv(AtmosModel & atms, const AerosolConcentration & aerocon
 	    double sha = 1 - (taer55p / aerocon.taer55);
 
 	    if (sha >= sham)
-		taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
+		taer55p = (double)(aerocon.taer55 * (1 - exp(-palt / 4)));
 	    else {
 		sha = -palt / log(sha);
-		taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / sha)));
+		taer55p = (double)(aerocon.taer55 * (1 - exp(-palt / sha)));
 	    }
 	}
     }

imagery/i.atcorr/altitude.h

@@ -41,41 +41,41 @@ struct AerosolConcentration;
 
 struct Altitude
 {
-	float xps;
-	float xpp;
+	double xps;
+	double xpp;
 
 	/* some vars */
-	mutable float palt;
-	float pps;
+	mutable double palt;
+	double pps;
 	int	  idatmp;
-	float taer55p;
-	float puw;
-	float puo3;
-	float ftray;
+	double taer55p;
+	double puw;
+	double puo3;
+	double ftray;
 
-	float puwus;
-	float puo3us;
+	double puwus;
+	double puo3us;
 
 	struct Plane_sim
 	{
-		float zpl[34];
-		float ppl[34];
-		float tpl[34];
-		float whpl[34];
-		float wopl[34];
+		double zpl[34];
+		double ppl[34];
+		double tpl[34];
+		double whpl[34];
+		double wopl[34];
 	} plane_sim;
 
 private:
     /* remember the original input values
      these values are set the first time when parse is called
      and used in subsequent calls to init to set xps and xpp */
-    float original_xps;
-    float original_xpp;
-    float original_taer55p;
-    float original_puw;
-    float original_puo3;
+    double original_xps;
+    double original_xpp;
+    double original_taer55p;
+    double original_puw;
+    double original_puo3;
 
-	void pressure(AtmosModel& atms, float& uw, float& uo3);
+	void pressure(AtmosModel& atms, double& uw, double& uo3);
 
 	void presplane(AtmosModel& atms);
 
@@ -86,7 +86,7 @@ public:
 	void print();
 
     /* Set the height to be used the next time init is called */
-    void set_height(const float height) { original_xps = height; }
+    void set_height(const double height) { original_xps = height; }
     /* call init only once: init parses input file */
     void init(AtmosModel& atms, const AerosolConcentration &aerocon);
     /* call update_hv whenever xps changes */

imagery/i.atcorr/atmosmodel.cpp

@@ -8,14 +8,14 @@ extern "C" {
 
 void AtmosModel::tropic()
 {
-    static const float z1[34] =
+    static const double 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 double 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,
@@ -23,7 +23,7 @@ void AtmosModel::tropic()
 	    3.05f, 1.59f, .854f, .0579f, 3e-4f, 0.f
 	};
 
-    static const float t1[34] =
+    static const double 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,
@@ -31,7 +31,7 @@ void AtmosModel::tropic()
 	    243.f, 254.f, 265.f, 270.f, 219.f, 210.f, 210.f
 	};
 
-    static const float wh1[34] =
+    static const double wh1[34] =
 	{ 
 	    19.f, 13.f, 9.3f, 4.7f, 2.2f, 1.5f, .85f, .47f, .25f, .12f, .05f, 
 
@@ -41,7 +41,7 @@ void AtmosModel::tropic()
 	    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 double 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,
@@ -63,14 +63,14 @@ void AtmosModel::tropic()
 
 void AtmosModel::midsum()
 {
-    static const float z1[34] =
+    static const double 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 double 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,
@@ -78,7 +78,7 @@ void AtmosModel::midsum()
 	    3.33f, 1.76f, .951f, .0671f, 3e-4f, 0.f
 	};
 
-    static const float t1[34] =
+    static const double 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,
@@ -86,7 +86,7 @@ void AtmosModel::midsum()
 	    270.f, 276.f, 218.f, 210.f, 210.f
 	};
 
-    static const float wh1[34] =
+    static const double 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,
@@ -94,7 +94,7 @@ void AtmosModel::midsum()
 	    1.1e-4f, 4.3e-5f, 1.9e-5f, 1.3e-6f, 1.4e-7f, 1e-9f, 0.f
 	};
 
-    static const float wo1[34] =
+    static const double 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,
@@ -116,14 +116,14 @@ void AtmosModel::midsum()
 
 void AtmosModel::midwin()
 {
-    static const float z1[34] =
+    static const double 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 double 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,
@@ -131,7 +131,7 @@ void AtmosModel::midwin()
 	    11.1f, 5.18f, 2.53f, 1.29f, .682f, .0467f, 3e-4f, 0.f
 	};
 
-    static const float t1[34] =
+    static const double 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,
@@ -139,7 +139,7 @@ void AtmosModel::midwin()
 	    215.2f, 217.4f, 227.8f, 243.2f, 258.5f, 265.7f, 230.7f, 210.2f, 210.f
 	};	
 
-    static const float wh1[34] =
+    static const double 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,
@@ -147,7 +147,7 @@ void AtmosModel::midwin()
 	    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 double wo1[34] = 
 	{ 
 	    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,
@@ -169,14 +169,14 @@ void AtmosModel::midwin()
 
 void AtmosModel::subsum()
 {
-    static const float z1[34] =
+    static const double 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 double 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,
@@ -184,7 +184,7 @@ void AtmosModel::subsum()
 	    3.4f, 1.81f, .987f, .0707f, 3e-4f, 0.f
 	};
 
-    static const float t1[34] =
+    static const double 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,
@@ -192,7 +192,7 @@ void AtmosModel::subsum()
 	    274.f, 277.f, 216.f, 210.f, 210.f
 	};
 
-    static const float wh1[34] =
+    static const double 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,
@@ -200,7 +200,7 @@ void AtmosModel::subsum()
 	    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 double wo1[34] = 
 	{ 
 	    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,
@@ -222,14 +222,14 @@ void AtmosModel::subsum()
 
 void AtmosModel::subwin()
 {
-    static const float z1[34] =
+    static const double 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 double 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,
@@ -237,7 +237,7 @@ void AtmosModel::subwin()
 	    22.56f, 10.2f, 4.701f, 2.243f, 1.113f, .5719f, .04016f, 3e-4f, 0.f
 	};
 
-    static const float t1[34] =
+    static const double 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,
@@ -245,7 +245,7 @@ void AtmosModel::subwin()
 	    211.2f, 216.f, 222.2f, 234.7f, 247.f, 259.3f, 245.7f, 210.f, 210.f
 	};
 
-    static const float wh1[34] =
+    static const double 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,
@@ -253,7 +253,7 @@ void AtmosModel::subwin()
 	    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 double 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,
@@ -275,14 +275,14 @@ void AtmosModel::subwin()
 
 void AtmosModel::us62()
 {
-    static const float z1[34] =
+    static const double 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 double 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,
@@ -290,7 +290,7 @@ void AtmosModel::us62()
 	    11.97f, 5.746f, 2.871f, 1.491f, .7978f, .0552f, 3.008e-4f, 0.f
 	};
 
-    static const float t1[34] =
+    static const double 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,
@@ -298,7 +298,7 @@ void AtmosModel::us62()
 	    221.6f, 226.5f, 236.5f, 253.4f, 264.2f, 270.6f, 219.7f, 210.f, 210.f
 	};
 
-    static const float wh1[34] =
+    static const double 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,
@@ -306,7 +306,7 @@ void AtmosModel::us62()
 	    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 double 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,

imagery/i.atcorr/atmosmodel.h

@@ -39,15 +39,15 @@ struct AtmosModel
 	long int idatm;	/* atmospheric model*/
 
 	/* secondary */
-    float uw;
-    float uo3;
+    double uw;
+    double uo3;
 
 	/* primary */
-	float z[34];
-	float p[34];
-	float t[34];
-	float wh[34];
-	float wo[34];
+	double z[34];
+	double p[34];
+	double t[34];
+	double wh[34];
+	double wo[34];
 
 private:
 	/* methods to initialize each model */

imagery/i.atcorr/common.cpp

@@ -1,6 +1,6 @@
 #include "common.h"
 #ifdef WIN32
-#pragma warning(disable:4305)	/* disable warning about initialization of a float by a double */
+#pragma warning(disable:4305)	/* disable warning about initialization of a double by a double */
 #endif
 
 Sixs_aer sixs_aer;	/* will be initialized in AerosolModel */

imagery/i.atcorr/common.h

@@ -38,60 +38,60 @@ using std::numeric_limits;
 const long int nt	= 26;
 
 /* Constants */
-const float sigma	= 0.056032f;
-const float delta	= 0.0279f;
-const float xacc	= 1.e-06f;
-const float step	= 0.0025f;
+const double sigma	= 0.056032f;
+const double delta	= 0.0279f;
+const double xacc	= 1.e-06f;
+const double step	= 0.0025f;
 
 
 /* Globals */
 /* not sure what the name stands for */
 struct Sixs_sos
 {
-	float phasel[10][83];
-	float cgaus[83];
-	float pdgs[83];
+	double phasel[10][83];
+	double cgaus[83];
+	double pdgs[83];
 };
 
 struct Sixs_aer
 {
-	float ext[10];
-	float ome[10]; 
-	float gasym[10]; 
-	float phase[10];
+	double ext[10];
+	double ome[10]; 
+	double gasym[10]; 
+	double phase[10];
 };
 
 struct Sixs_aerbas
 {
-	float bdm_ph[10][83];		/* background desert model... */
-	float bbm_ph[10][83];		/* biomass burning model... */
-	float stm_ph[10][83];		/* stratospherique aerosol model... */
-	float dust_ph[10][83];		/* dust model */
-	float wate_ph[10][83];		/* water model */
-	float ocea_ph[10][83];		/* ocean model */
-	float soot_ph[10][83];		/* soot model */
-
-	float usr_ph[10][83];		/* user defined model from size distribution */
-	float (*ph)[10][83];		/* pointer to current active model */
+	double bdm_ph[10][83];		/* background desert model... */
+	double bbm_ph[10][83];		/* biomass burning model... */
+	double stm_ph[10][83];		/* stratospherique aerosol model... */
+	double dust_ph[10][83];		/* dust model */
+	double wate_ph[10][83];		/* water model */
+	double ocea_ph[10][83];		/* ocean model */
+	double soot_ph[10][83];		/* soot model */
+
+	double usr_ph[10][83];		/* user defined model from size distribution */
+	double (*ph)[10][83];		/* pointer to current active model */
 };
 
 struct Sixs_trunc
 {
-	float pha[83];
-	float betal[81];
+	double pha[83];
+	double betal[81];
 };
 
 struct Sixs_disc
 {
-	float roatm[3][10];
-	float dtdir[3][10];
-	float dtdif[3][10];
-	float utdir[3][10];
-	float utdif[3][10];
-	float sphal[3][10];
-	float wldis[10];
-	float trayl[10];
-    float traypl[10];
+	double roatm[3][10];
+	double dtdir[3][10];
+	double dtdif[3][10];
+	double utdir[3][10];
+	double utdif[3][10];
+	double sphal[3][10];
+	double wldis[10];
+	double trayl[10];
+    double traypl[10];
 };
 
 extern Sixs_sos sixs_sos;

imagery/i.atcorr/gauss.cpp

@@ -4,8 +4,8 @@
 
 const long int mu2	= 48;
 
-float Gauss::angmu[10] = { 85.f, 80.f, 70.f, 60.f, 50.f, 40.f, 30.f, 20.f, 10.f, 0.f };
-float Gauss::angphi[13] = { 0.f, 30.f, 60.f, 90.f, 120.f, 150.f, 180.f, 210.f, 240.f, 270.f, 300.f, 330.f, 360.f };
+double Gauss::angmu[10] = { 85.f, 80.f, 70.f, 60.f, 50.f, 40.f, 30.f, 20.f, 10.f, 0.f };
+double Gauss::angphi[13] = { 0.f, 30.f, 60.f, 90.f, 120.f, 150.f, 180.f, 210.f, 240.f, 270.f, 300.f, 330.f, 360.f };
 
 /*  preliminary computations for gauss integration */
 void Gauss::init()
@@ -13,13 +13,13 @@ void Gauss::init()
     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);
+    for (j = 0; j < 13; ++j) angphi[j] = (double)(angphi[j] * M_PI / 180.f);
+    for (j = 0; j < 10; ++j) angmu[j] =	 (double)cos(angmu[j] * M_PI / 180.f);
 
     /* calculate rm & gb */
 	
-    float anglem[mu2];
-    float weightm[mu2];
+    double anglem[mu2];
+    double weightm[mu2];
     gauss (-1.f, 1.f, anglem, weightm, mu2);
 
     gb[STDI(-mu)]   = 0;
@@ -42,14 +42,14 @@ void Gauss::init()
     }
 
     /* calculate rp & gp */
-    gauss (0.f, (float)2 * M_PI, rp, gp, np);
+    gauss (0.f, (double)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. */
-void Gauss::gauss (float a, float b, float *x, float *w, long int n)
+void Gauss::gauss (double a, double b, double *x, double *w, long int n)
 {
     int m = (n + 1) / 2;
     double xm = (b + a) / 2;
@@ -79,9 +79,9 @@ void Gauss::gauss (float a, float b, float *x, float *w, long int n)
 	} while(fabs(z - z1) > 3e-14);
 
 	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));
+	x[i] = (double)(xm - xl * z);
+        x[n - 1 - i] = (double)(xm + xl * z);
+	w[i] = (double)(2 * xl / ((1 - z * z) * pp * pp));
 	w[n - 1 - i] = w[i];
     }
 }

imagery/i.atcorr/gauss.h

@@ -10,17 +10,17 @@ const long int np	= 49;
 struct Gauss
 {
 private:
-	static float angmu[10];
-	static float angphi[13];
+	static double angmu[10];
+	static double angphi[13];
 
 public:
 	/* [a,b] = [0,2*Pi] */
-	float rp[np];			/* gaussian angles */
-	float gp[np];			/* gaussian weights */
+	double rp[np];			/* gaussian angles */
+	double gp[np];			/* gaussian weights */
 
 	// [a,b] = [-1,1]
-	float rm[2*mu+1];		/* shifted gaussian angles */
-	float gb[2*mu+1];		/* shifted gaussian weights */
+	double rm[2*mu+1];		/* shifted gaussian angles */
+	double gb[2*mu+1];		/* shifted gaussian weights */
 					/* with the ends zeroed as well as the center */
 					/* [0 ? ? ? ? 0 ? ? ? ? 0] */
 
@@ -30,7 +30,7 @@ public:
 	/*	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. */
-	static void gauss (float a, float b, float *x, float *w, long int n);
+	static void gauss (double a, double b, double *x, double *w, long int n);
 };
 
 #endif /* MY_GAUSS_H */

imagery/i.atcorr/geomcond.cpp

@@ -46,7 +46,7 @@ extern "C" {
   return dsol		
   dsol is a multiplicative factor to apply to the mean value of solar constant 
 */
-float GeomCond::varsol ()
+double GeomCond::varsol ()
 {
 /* calculation of the variability of the solar constant during the year. 
    jday is the number of the day in the month   */
@@ -56,13 +56,13 @@ float GeomCond::varsol ()
     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 ((double) (j - 4) * 0.9856f * M_PI / 180.f) * .01673f;
+    return 1.f / (double)(tmp * tmp);
 }
 
 
 /* spot, landsat5 and landsat7 is handled the same way */
-void GeomCond::landsat(float tu)
+void GeomCond::landsat(double tu)
 {
 /*     warning !!! */
 /*     xlon and xlat are the coordinates of the scene center. */
@@ -77,7 +77,7 @@ void GeomCond::landsat(float tu)
   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)
+void GeomCond::possol(double tu)
 {
     long int ia = 0;
     long int nojour;
@@ -107,7 +107,7 @@ void GeomCond::day_number(long int ia, long int& j)
 /* returns the sign of the element */
 #define SIGN(X) (((X) >= 0) ? 1. : -1.) 
 
-void GeomCond::pos_fft (long int j, float tu)
+void GeomCond::pos_fft (long int j, double tu)
 {
     /* Local variables */
     double ah, et, az, caz, xla, tet, tsm, tsv, elev, azim, gdelta, amuzero;
@@ -119,7 +119,7 @@ void GeomCond::pos_fft (long int j, float tu)
     /* mean solar time (heure decimale) */
     tsm = tu + xlon / 15.;
     xla = xlat * M_PI / 180.;
-    tet = (float)(j) * M_PI2 / 365.;
+    tet = (double)(j) * M_PI2 / 365.;
 
     /* time equation (in mn.dec) */
     et = 7.5e-5f + 0.001868f * cos (tet) - 0.032077f * sin (tet) - 
@@ -158,8 +158,8 @@ void GeomCond::pos_fft (long int j, float tu)
     elev = elev * 180. / M_PI;
 	
     /*     conversion in degrees */
-    asol = (float)(90. - elev);
-    phi0 = (float)(azim * 180. / M_PI);
+    asol = (double)(90. - elev);
+    phi0 = (double)(azim * 180. / M_PI);
 }
 
 /*
@@ -168,7 +168,7 @@ void GeomCond::pos_fft (long int j, float tu)
   2 = goes east observation
   3 = goes west observation
 */
-void GeomCond::posobs(float tu, int nc, int nl)
+void GeomCond::posobs(double tu, int nc, int nl)
 {
     double yr, xr, alti;
 
@@ -239,11 +239,11 @@ void GeomCond::posobs(float tu, int nc, int nl)
 	ylon = atan(xt / zt);
     }
  
-    xlat = (float)(ylat * crd);
+    xlat = (double)(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 = (double)(ylon * crd);
+    else if(igeom == 2) xlon = (double)(ylon * crd - 75.);
+    else xlon = (double)(ylon * crd - 135.);
  
     possol(tu);
  
@@ -253,11 +253,11 @@ void GeomCond::posobs(float tu, int nc, int nl)
 
     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);
+    avis = (double)(asin((1. + alti / re) * (gam)) * 180. / M_PI);
+    phiv = (double)((atan2(tan(ylon),sin(ylat)) + M_PI) * 180. / M_PI);
 }
 
-void GeomCond::posnoa(float tu, int nc, float xlonan, float campm, float hna)
+void GeomCond::posnoa(double tu, int nc, double xlonan, double campm, double hna)
 {
 /*     noaa 6 definition
        orbite inclination ai in radians
@@ -276,7 +276,7 @@ void GeomCond::posnoa(float tu, int nc, float xlonan, float campm, float hna)
     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));
+    avis = (double)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);
@@ -291,8 +291,8 @@ void GeomCond::posnoa(float tu, int nc, float xlonan, float campm, float hna)
 	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);
+    xlat = (double)(ylat * 180. / M_PI);
+    xlon = (double)(ylo1 * 180. / M_PI);
  
 
 
@@ -304,11 +304,11 @@ void GeomCond::posnoa(float tu, int nc, float xlonan, float campm, float hna)
     {
 	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);
+	phiv = (double)atan2(xnum,xden);
     }
     else phiv = 0.;
-    phiv = (float)(phiv * 180. / M_PI);
-    avis = (float)(fabs(avis) * 180. / M_PI);
+    phiv = (double)(phiv * 180. / M_PI);
+    avis = (double)(fabs(avis) * 180. / M_PI);
 }
 
 void GeomCond::parse()
@@ -316,8 +316,8 @@ void GeomCond::parse()
     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;
+    double campm = -1.0f;	/* initialize in case igeom == 5 */
+    double tu, xlonan, hna;
     int nc, nl;
 
     switch(igeom)
@@ -346,7 +346,7 @@ void GeomCond::parse()
 	posobs(tu, nc, nl);
 	break;
     }
-    case 4: campm = 1.0f; /* avhrr PM */
+    case 4: campm = 1.0f; /* avhrr PM, case 4 must fall through to case 5 */
     case 5:  		  /* avhrr PM and avhrr AM */
     {
 	cin >> month;
@@ -407,20 +407,20 @@ void GeomCond::parse()
     /*    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;
+    phi = (double)fabs(phiv - phi0);
+    phirad = (phi0 - phiv) * (double)M_PI / 180.f;
+    if (phirad < 0.f) phirad += (double)M_PI2;
+    if (phirad > M_PI2) phirad -= (double)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;
+    xmus = (double)cos (asol * M_PI / 180.f);
+    xmuv = (double)cos (avis * M_PI / 180.f);
+    xmup = (double)cos (phirad);
+    xmud = -xmus * xmuv - (double)sqrt (1.f - xmus * xmus) * (double)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;
+    adif = (double)acos (xmud) * 180.f / (double)M_PI;
 
     dsol = varsol();
 }

imagery/i.atcorr/geomcond.h

@@ -101,39 +101,39 @@ struct GeomCond
 	long int igeom;	/* geometrical conditions */
 
 	/* primary */
-	float asol;
-	float phi0;
-	float avis;
-	float phiv;
+	double asol;
+	double phi0;
+	double avis;
+	double phiv;
 	long int month;
 	long int jday;
-	float xlon;
-	float xlat;
+	double xlon;
+	double xlat;
 
 	/* some vars */
-	float phi;
-	float phirad;
-	float xmus; 
-	float xmuv; 
-	float xmup; 
-	float xmud;
-	float adif;
+	double phi;
+	double phirad;
+	double xmus; 
+	double xmuv; 
+	double xmup; 
+	double xmud;
+	double adif;
 
-    float dsol;
+    double dsol;
 
 	void  print();
 
 private:
 	/* conversion routines */
-	void possol(float tu);
-	void landsat(float tu);
-	void posobs(float tu, int nc, int nl);
-	void posnoa(float tu, int nc, float xlonan, float campm, float hna);
+	void possol(double tu);
+	void landsat(double tu);
+	void posobs(double tu, int nc, int nl);
+	void posnoa(double tu, int nc, double xlonan, double campm, double hna);
 
 	void day_number(long int ia, long int& j);
-	void pos_fft (long int j, float tu);
+	void pos_fft (long int j, double tu);
 
-	float varsol();	/* returns dsol as in fortran proggie */
+	double varsol();	/* returns dsol as in fortran proggie */
 	void parse();
 public:
 	static GeomCond Parse();

imagery/i.atcorr/interp.cpp

@@ -2,8 +2,8 @@
 #include "interp.h"
 
 void interp (const int iaer, const int idatmp, 
-	     const float wl, const float taer55, 
-	     const float taer55p, const float xmud, 
+	     const double wl, const double taer55, 
+	     const double taer55p, const double xmud, 
 	     InterpStruct& is)
 {
 /*     that for the atmosphere :
@@ -42,12 +42,12 @@ void interp (const int iaer, const int idatmp,
     /*    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;
+    double alphaa = 0;
+    double betaa = 0;
+    double alphar = 0;
+    double betar = 0;
+    double alphac = 0;
+    double betac = 0;
     is.phaa = 0;
     is.roaero = 0;
     is.dtota = 1;
@@ -55,17 +55,17 @@ void interp (const int iaer, const int idatmp,
     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];
+    double coef = (double)log(sixs_disc.wldis[lsup] / sixs_disc.wldis[linf]);
+    double 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 = (double)(log(sixs_aer.phase[lsup] / sixs_aer.phase[linf]) / coef);
+	betaa = (double)(sixs_aer.phase[linf] / pow(wlinf,alphaa));
+	is.phaa = (double)(betaa * pow(wl,alphaa));
     }
 
-    float d2 = 2 + delta;
+    double d2 = 2 + delta;
     is.phar = (2 * (1 - delta) / d2) * .75f * (1 + xmud * xmud) + 3 * delta / d2;
     if(idatmp == 0)
     {
@@ -82,9 +82,9 @@ void interp (const int iaer, const int idatmp,
 	}
 	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));
+	    alphar = (double)(log(sixs_disc.roatm[0][lsup] / sixs_disc.roatm[0][linf]) / coef);
+	    betar = (double)(sixs_disc.roatm[0][linf] / pow(wlinf,alphar));
+	    is.rorayl = (double)(betar * pow(wl,alphar));
 	}
 
 	if(sixs_disc.roatm[1][linf] < 0.001)
@@ -94,9 +94,9 @@ void interp (const int iaer, const int idatmp,
 	}
         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));
+	    alphac = (double)(log(sixs_disc.roatm[1][lsup] / sixs_disc.roatm[1][linf]) / coef);
+	    betac = (double)(sixs_disc.roatm[1][linf] / pow(wlinf,alphac));
+	    is.romix = (double)(betac * pow(wl,alphac));
 	}
 
 	if(iaer != 0)
@@ -109,95 +109,95 @@ void interp (const int iaer, const int idatmp,
 	    }
 	    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));
+		alphaa = (double)(log(sixs_disc.roatm[2][lsup] / sixs_disc.roatm[2][linf]) / coef);
+		betaa = (double)(sixs_disc.roatm[2][linf] / pow(wlinf,alphaa));
+		is.roaero = (double)(betaa * pow(wl,alphaa));
 	    }
 	}
     }
 
-    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));
+    alphar = (double)(log(sixs_disc.trayl[lsup] / sixs_disc.trayl[linf]) / coef);
+    betar = (double)(sixs_disc.trayl[linf] / pow(wlinf,alphar));
+    is.tray = (double)(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));
+	alphar = (double)(log(sixs_disc.traypl[lsup] / sixs_disc.traypl[linf]) / coef);
+        betar = (double)(sixs_disc.traypl[linf] / pow(wlinf,alphar));
+        is.trayp = (double)(betar * pow(wl,alphar));
     }
     else is.trayp = 0;
 
     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]);
+	alphaa = (double)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] / (sixs_aer.ext[linf] * sixs_aer.ome[linf])) / coef);
+	betaa = (double)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
+	is.tsca = (double)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+	alphaa = (double)(log(sixs_aer.ext[lsup] / sixs_aer.ext[linf]) / coef);
+	betaa = (double)(sixs_aer.ext[linf] / pow(wlinf,alphaa));
+	is.taerp = (double)(taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+	is.taer = (double)(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];
+    double drinf = sixs_disc.dtdif[0][linf] + sixs_disc.dtdir[0][linf];
+    double drsup = sixs_disc.dtdif[0][lsup] + sixs_disc.dtdir[0][lsup];
+    alphar = (double)(log(drsup / drinf) / coef);
+    betar = (double)(drinf / pow(wlinf,alphar));
+    is.dtotr = (double)(betar * pow(wl,alphar));
+    double dtinf = sixs_disc.dtdif[1][linf] + sixs_disc.dtdir[1][linf];
+    double dtsup = sixs_disc.dtdif[1][lsup] + sixs_disc.dtdir[1][lsup];
+    alphac = (double)(log((dtsup * drinf) / (dtinf * drsup)) / coef);
+    betac = (double)((dtinf / drinf) / pow(wlinf,alphac));
+    double dtotc = (double)(betac * pow(wl,alphac));
+    double dainf = sixs_disc.dtdif[2][linf] + sixs_disc.dtdir[2][linf];
+    double 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));
+	alphaa = (double)(log(dasup / dainf) / coef);
+	betaa = (double)(dainf / pow(wlinf,alphaa));
+	is.dtota = (double)(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];
+    double urinf = sixs_disc.utdif[0][linf] + sixs_disc.utdir[0][linf];
+    double ursup = sixs_disc.utdif[0][lsup] + sixs_disc.utdir[0][lsup];
+    alphar = (double)(log(ursup / urinf) / coef);
+    betar = (double)(urinf / pow(wlinf,alphar));
+    is.utotr = (double)(betar * pow(wl,alphar));
+    double utinf = sixs_disc.utdif[1][linf] + sixs_disc.utdir[1][linf];
+    double utsup = sixs_disc.utdif[1][lsup] + sixs_disc.utdir[1][lsup];
+    alphac = (double)(log((utsup * urinf) / (utinf * ursup)) / coef);
+    betac = (double)((utinf / urinf) / pow(wlinf,alphac));
+    double utotc = (double)(betac * pow(wl,alphac));
+    double uainf = sixs_disc.utdif[2][linf] + sixs_disc.utdir[2][linf];
+    double uasup = sixs_disc.utdif[2][lsup] + sixs_disc.utdir[2][lsup];
     is.utott = utotc * is.utotr;
 
     if(iaer != 0)
     {
-	alphaa = (float)(log(uasup / uainf) / coef);
-	betaa = (float)(uainf / pow(wlinf,alphaa));
-	is.utota = (float)(betaa * pow(wl,alphaa));
+	alphaa = (double)(log(uasup / uainf) / coef);
+	betaa = (double)(uainf / pow(wlinf,alphaa));
+	is.utota = (double)(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];
+    double arinf = sixs_disc.sphal[0][linf];
+    double arsup = sixs_disc.sphal[0][lsup];
+    alphar = (double)(log(arsup / arinf) / coef);
+    betar = (double)(arinf / pow(wlinf,alphar));
+    is.asray = (double)(betar * pow(wl,alphar));
+    double atinf = sixs_disc.sphal[1][linf];
+    double atsup = sixs_disc.sphal[1][lsup];
+    alphac = (double)(log(atsup / atinf) / coef);
+    betac = (double)(atinf / pow(wlinf,alphac));
+    is.astot = (double)(betac * pow(wl,alphac));
+    double aainf = sixs_disc.sphal[2][linf];
+    double 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));
+	alphaa = (double)(log(aasup / aainf) / coef);
+	betaa = (double)(aainf / pow(wlinf,alphaa));
+	is.asaer = (double)(betaa * pow(wl,alphaa));
     }
 }

imagery/i.atcorr/interp.h

@@ -3,25 +3,25 @@
 
 struct InterpStruct
 {
-	float romix; 
-	float rorayl; 
-	float roaero;
-	float phaa; 
-	float phar; 
-	float tsca;
-	float tray; 
-	float trayp; 
-	float taer;
-	float taerp; 
-	float dtott; 
-	float utott;
-	float astot; 
-	float asray; 
-	float asaer;
-	float utotr; 
-	float utota; 
-	float dtotr;
-	float dtota;
+	double romix; 
+	double rorayl; 
+	double roaero;
+	double phaa; 
+	double phar; 
+	double tsca;
+	double tray; 
+	double trayp; 
+	double taer;
+	double taerp; 
+	double dtott; 
+	double utott;
+	double astot; 
+	double asray; 
+	double asaer;
+	double utotr; 
+	double utota; 
+	double dtotr;
+	double dtota;
 };
 
 /*
@@ -29,8 +29,8 @@ To estimate the different atmospheric functions r(mS,mv,fS,fv), T(q) and S at an
 wavelength from the 10 discret computations (subroutine DISCOM).
  */
 void interp (const int iaer, const int idatmp, 
-			 const float wl, const float taer55, 
-			 const float taer55p, const float xmud, 
+			 const double wl, const double taer55, 
+			 const double taer55p, const double xmud, 
 			 InterpStruct& is);
 
 #endif /* INTERP_H */

imagery/i.atcorr/iwave.cpp

@@ -1463,7 +1463,7 @@ void IWave::etmplus(int iwa)
     }
 }
 
-float IWave::solirr(const float fwl) const
+double IWave::solirr(const double fwl) const
 {
 /*    si (in w/m2/micron) contains the values of the solar
       irradiance between 0.25 and 4.0 microns, by step of 0.0025 m.
@@ -1688,7 +1688,7 @@ float IWave::solirr(const float fwl) const
 	8.85,   8.83,   8.81
     };
 
-    float pas = 0.0025;
+    double pas = 0.0025;
     int   iwl = (int)((fwl - 0.250) / pas + 1.5);
 	  
     if(iwl >= 0) return si[iwl-1];
@@ -4862,18 +4862,18 @@ void IWave::planetscope0f10(int iwa)
 
 /* filter functions must be defined above */
 
-float IWave::equivwl() const
+double IWave::equivwl() const
 {
-    float seb = 0;
-    float wlwave = 0;
+    double seb = 0;
+    double wlwave = 0;
 
     for(int i = iinf; i <= isup; i++)
     {
-	float sbor = ffu.s[i];
+	double sbor = ffu.s[i];
 	if(i == iinf || i == isup) sbor *= 0.5;
-	float fwl = (float)(0.25 + i * step);
-	float swl = solirr(fwl);
-	float coef = sbor * step * swl;
+	double fwl = (double)(0.25 + i * step);
+	double swl = solirr(fwl);
+	double coef = sbor * step * swl;
 	seb += coef;
 	wlwave += fwl * coef;
     }
@@ -4952,7 +4952,7 @@ void IWave::parse()
 
 	if (iwave > 1) {
 	    int imax;
-	    float smax, sthreshold;
+	    double smax, sthreshold;
 
 	    imax = -1;
 	    smax = 0;

imagery/i.atcorr/iwave.h

@@ -42,15 +42,15 @@ struct IWave
 	int iinf;
 	int isup;
 
-	float wl;
-	float wlmoy;
+	double wl;
+	double wlmoy;
 
 	
 	struct FFu
 	{
-		float s[1501];
-		float wlinf;
-		float wlsup;
+		double s[1501];
+		double wlinf;
+		double wlsup;
 	} ffu;
 
 private:	
@@ -93,12 +93,12 @@ public:
 	/* To compute the equivalent wavelength needed for the calculation of the
 	  downward radiation field used in the computation of the non lambertian 
 	  target contribution (main.f). */
-	float equivwl() const;
+	double equivwl() const;
 
 	/* To read the solar irradiance (in Wm-2mm-1) from 250 nm to 4000 nm by 
 	steps of 2.5 nm, The total solar irradiance is put equal to 1372 Wm-2. 
 	Between 250 and 4000 nm we have 1358 Wm-2. */
-	float solirr(float wl) const;
+	double solirr(double wl) const;
 
 	void print();
 	static IWave Parse();

imagery/i.atcorr/main.cpp

@@ -138,7 +138,7 @@ static CELL round_c(FCELL x)
 }
 
 /* round height, input and output unit is meters */
-static int round_h(float x)
+static int round_h(double x)
 {
     x /= BIN_ALT;
     x = floor(x + 0.5);
@@ -148,7 +148,7 @@ static int round_h(float x)
 }
 
 /* round visibility, input unit is km, output unit is meters */
-static int round_v(float x)
+static int round_v(double x)
 {
     x *= 1000. / BIN_VIS;
     x = floor(x + 0.5);
@@ -205,7 +205,7 @@ class TICache
     unsigned int tree_size;
 
   private:
-    struct RBitem set_alt_vis(float alt, float vis)
+    struct RBitem set_alt_vis(double alt, double vis)
     {
 	struct RBitem rbitem;
 
@@ -222,7 +222,7 @@ class TICache
 	RBTree = rbtree_create(compare_hv, sizeof(struct RBitem));
 	tree_size = 0;
     }
-    int search(float alt, float vis, TransformInput *ti)
+    int search(double alt, double vis, TransformInput *ti)
     {
 	struct RBitem search_ti = set_alt_vis(alt, vis);
 	struct RBitem *found_ti;
@@ -237,7 +237,7 @@ class TICache
 
     }
 
-    void add(TransformInput ti, float alt, float vis)
+    void add(TransformInput ti, double alt, double vis)
     {
 	struct RBitem insert_ti = set_alt_vis(alt, vis);
 
@@ -402,15 +402,17 @@ static void process_raster(int ifd, InputMask imask, ScaleRange iscale,
 	    G_debug(3, "Computed r%d (%d), c%d (%d)", row, nrows, col, ncols);
 	    /* transform from iscale.[min,max] to [0,1] */
 	    buf[col] =
-		(buf[col] - iscale.min) / ((float)iscale.max -
-					   (float)iscale.min);
+		(buf[col] - iscale.min) / ((double)iscale.max -
+					   (double)iscale.min);
 	    buf[col] = transform(ti, imask, buf[col]);
+	    if (Rast_is_f_null_value(&buf[col]))
+		G_fatal_error(_("Numerical instability in 6S"));
 	    /* transform from [0,1] to oscale.[min,max] */
 	    buf[col] =
-		buf[col] * ((float)oscale.max - (float)oscale.min) +
+		buf[col] * ((double)oscale.max - (double)oscale.min) +
 		oscale.min;
 
-	    if (oint && (buf[col] > (float)oscale.max))
+	    if (oint && (buf[col] > (double)oscale.max))
 		G_warning(_("The output data will overflow. Reflectance > 100%%"));
 	}
 
@@ -619,7 +621,7 @@ int main(int argc, char *argv[])
 			  opts.ivis->answer);
     }
 
-    /* open a floating point raster or not? */
+    /* open a doubleing point raster or not? */
     if (opts.oint->answer) {
 	if ((oimg_fd = Rast_open_new(opts.oimg->answer, CELL_TYPE)) < 0)
 	    G_fatal_error(_("Unable to create raster map <%s>"),

imagery/i.atcorr/transform.cpp

@@ -1,7 +1,7 @@
 #include <math.h>
 #include "transform.h"
 
-void EtmDN(int iwave, float asol, bool before, float &lmin, float &lmax)
+void EtmDN(int iwave, double asol, bool before, double &lmin, double &lmax)
 {
     if (before)		/* ETM+ digital numbers taken before July 1, 2000 */
     {
@@ -132,13 +132,13 @@ void EtmDN(int iwave, float asol, bool before, float &lmin, float &lmax)
 /* 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 */
-float transform(const TransformInput ti, InputMask imask, float idn)
+double transform(const TransformInput ti, InputMask imask, double idn)
 {
     /* convert from radiance to reflectance */
     if((imask & ETM_BEFORE) || (imask & ETM_AFTER))
     {
         /* http://ltpwww.gsfc.nas */
-        float lmin, lmax;
+        double lmin, lmax;
         EtmDN(ti.iwave, ti.asol, imask & ETM_BEFORE, lmin, lmax);
 
         /* multiply idn by 255.f to correct precondition that idn lies in [0, 255] */
@@ -146,16 +146,16 @@ float transform(const TransformInput ti, InputMask imask, float idn)
         if (idn < 0.f) idn = 0.f;
         idn /= 255.f;
     }
-    if(imask & RADIANCE) idn += (float)M_PI * idn * 255.f * ti.sb / ti.xmus / ti.seb;
+    if(imask & RADIANCE) idn += (double)M_PI * idn * 255.f * ti.sb / ti.xmus / ti.seb;
           
-    float rapp = idn;
-    float ainrpix = ti.ainr[0][0];
+    double rapp = idn;
+    double ainrpix = ti.ainr[0][0];
     /*
-    float xa = 0.0f;
-    float xb = 0.0f;
-    float xc = 0.0f;
+    double xa = 0.0f;
+    double xb = 0.0f;
+    double xc = 0.0f;
     */
-    float rog = rapp / ti.tgasm;
+    double 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
@@ -165,7 +165,7 @@ float transform(const TransformInput ti, InputMask imask, float idn)
        bright pixels in the image. Check the output file to see if that
        has happened. */
 
-    float decrfact = 1.0f;
+    double decrfact = 1.0f;
     if (rog < (ainrpix / ti.tgasm))
     {
 	do
@@ -179,7 +179,7 @@ float transform(const TransformInput ti, InputMask imask, float idn)
     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;
+    xa = (double)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;
     */

imagery/i.atcorr/transform.h

@@ -13,17 +13,17 @@
 struct TransformInput
 {
     int iwave;
-    float asol;
+    double asol;
     
-    float ainr[2][3];
-    float sb;
-    float seb;
-    float tgasm;
-    float sutott;
-    float sdtott;
-    float sast;
-    float srotot;
-    float xmus;
+    double ainr[2][3];
+    double sb;
+    double seb;
+    double tgasm;
+    double sutott;
+    double sdtott;
+    double sast;
+    double srotot;
+    double xmus;
 };
 
 /* The following combinations of input values types exist */
@@ -42,6 +42,6 @@ enum InputMask
 /* 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 */
-extern float transform(const TransformInput ti, InputMask imask, float idn);
+extern double transform(const TransformInput ti, InputMask imask, double idn);
 
 #endif /* TRANSFORM_H */