[GRASS-SVN] r73183 - grass/trunk/imagery/i.atcorr
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Aug 25 14:42:22 PDT 2018
Author: mmetz
Date: 2018-08-25 14:42:21 -0700 (Sat, 25 Aug 2018)
New Revision: 73183
Modified:
grass/trunk/imagery/i.atcorr/6s.cpp
grass/trunk/imagery/i.atcorr/6s.h
grass/trunk/imagery/i.atcorr/abstra.cpp
grass/trunk/imagery/i.atcorr/abstra.h
grass/trunk/imagery/i.atcorr/aerosolconcentration.cpp
grass/trunk/imagery/i.atcorr/aerosolconcentration.h
grass/trunk/imagery/i.atcorr/aerosolmodel.cpp
grass/trunk/imagery/i.atcorr/aerosolmodel.h
grass/trunk/imagery/i.atcorr/altitude.cpp
grass/trunk/imagery/i.atcorr/altitude.h
grass/trunk/imagery/i.atcorr/atmosmodel.cpp
grass/trunk/imagery/i.atcorr/atmosmodel.h
grass/trunk/imagery/i.atcorr/common.cpp
grass/trunk/imagery/i.atcorr/common.h
grass/trunk/imagery/i.atcorr/computations.cpp
grass/trunk/imagery/i.atcorr/gauss.cpp
grass/trunk/imagery/i.atcorr/gauss.h
grass/trunk/imagery/i.atcorr/geomcond.cpp
grass/trunk/imagery/i.atcorr/geomcond.h
grass/trunk/imagery/i.atcorr/interp.cpp
grass/trunk/imagery/i.atcorr/interp.h
grass/trunk/imagery/i.atcorr/iwave.cpp
grass/trunk/imagery/i.atcorr/iwave.h
grass/trunk/imagery/i.atcorr/main.cpp
grass/trunk/imagery/i.atcorr/transform.cpp
grass/trunk/imagery/i.atcorr/transform.h
Log:
i.atcorr: avoid float to double to float conversion, causing numerical instability
Modified: grass/trunk/imagery/i.atcorr/6s.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/6s.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/6s.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -22,10 +22,10 @@
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,7 +91,7 @@
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;
@@ -98,7 +98,7 @@
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,7 +107,7 @@
}
/* 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);
@@ -114,38 +114,38 @@
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 @@
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 @@
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;
+ 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;
- 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;
+ 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;
- float ani[2][3];
- float aini[2][3];
- float anr[2][3];
- float ainr[2][3];
+ 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 @@
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 @@
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 @@
}
- 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 @@
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 @@
}
/* 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 @@
srotot = srotot / seb;
alumet = alumet / sb;
/*
- float pizera = 0.0f;
+ double pizera = 0.0f;
if(aero.iaer != 0) pizera = ssdaer / sodaer;
*/
sodray = sodray / seb;
Modified: grass/trunk/imagery/i.atcorr/6s.h
===================================================================
--- grass/trunk/imagery/i.atcorr/6s.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/6s.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -29,9 +29,9 @@
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
Modified: grass/trunk/imagery/i.atcorr/abstra.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/abstra.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/abstra.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 */
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 */
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 */
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 */
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 */
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 */
-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 @@
} /* 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 @@
} /* 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
-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 @@
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 */
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 */
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 */
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 */
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 */
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 @@
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 */
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 */
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 */
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 */
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 */
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 @@
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 */
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 */
-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 @@
} /* 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 */
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 */
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 @@
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 @@
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 @@
} /* 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 */
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 @@
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 @@
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 @@
}
/* 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 @@
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 @@
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 @@
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 @@
{
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 @@
}
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 @@
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 @@
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 @@
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 @@
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];
Modified: grass/trunk/imagery/i.atcorr/abstra.h
===================================================================
--- grass/trunk/imagery/i.atcorr/abstra.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/abstra.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
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 */
Modified: grass/trunk/imagery/i.atcorr/aerosolconcentration.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/aerosolconcentration.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/aerosolconcentration.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -32,17 +32,17 @@
{
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 @@
/* 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 @@
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;
}
}
Modified: grass/trunk/imagery/i.atcorr/aerosolconcentration.h
===================================================================
--- grass/trunk/imagery/i.atcorr/aerosolconcentration.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/aerosolconcentration.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -32,17 +32,17 @@
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);
};
Modified: grass/trunk/imagery/i.atcorr/aerosolmodel.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/aerosolmodel.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/aerosolmodel.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -7,7 +7,7 @@
#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 @@
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 @@
{
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 @@
{
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 @@
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 @@
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 @@
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,31 +535,31 @@
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
@@ -566,13 +566,13 @@
};
- 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 @@
{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 @@
{
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 @@
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 @@
{
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 @@
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 @@
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 @@
}
}
-AerosolModel AerosolModel::Parse(const float xmud)
+AerosolModel AerosolModel::Parse(const double xmud)
{
AerosolModel aero;
aero.parse(xmud);
Modified: grass/trunk/imagery/i.atcorr/aerosolmodel.h
===================================================================
--- grass/trunk/imagery/i.atcorr/aerosolmodel.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/aerosolmodel.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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,16 +110,16 @@
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;
@@ -126,18 +126,18 @@
};
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);
};
Modified: grass/trunk/imagery/i.atcorr/altitude.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/altitude.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/altitude.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
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 @@
/* 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 @@
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 @@
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;
+ 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;
- float rmo3[34];
- float rmwh[34];
+ 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 @@
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 @@
xps = original_xps;
xpp = original_xpp;
- float uwus;
- float uo3us;
+ double uwus;
+ double uo3us;
if(xps <= 0)
{
xps = 0;
@@ -265,7 +265,7 @@
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 @@
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 @@
xps = original_xps;
xpp = original_xpp;
- float uwus;
- float uo3us;
+ double uwus;
+ double uo3us;
if (xps <= 0) {
xps = 0;
@@ -346,7 +346,7 @@
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 @@
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)));
}
}
}
Modified: grass/trunk/imagery/i.atcorr/altitude.h
===================================================================
--- grass/trunk/imagery/i.atcorr/altitude.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/altitude.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -41,28 +41,28 @@
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:
@@ -69,13 +69,13 @@
/* 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 @@
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 */
Modified: grass/trunk/imagery/i.atcorr/atmosmodel.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/atmosmodel.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/atmosmodel.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -8,7 +8,7 @@
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,
@@ -15,7 +15,7 @@
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 @@
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 @@
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 @@
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,7 +63,7 @@
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,
@@ -70,7 +70,7 @@
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 @@
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 @@
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 @@
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,7 +116,7 @@
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,
@@ -123,7 +123,7 @@
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 @@
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 @@
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 @@
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,7 +169,7 @@
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,
@@ -176,7 +176,7 @@
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 @@
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 @@
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 @@
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,7 +222,7 @@
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,
@@ -229,7 +229,7 @@
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 @@
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 @@
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 @@
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,7 +275,7 @@
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,
@@ -282,7 +282,7 @@
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 @@
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 @@
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 @@
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,
Modified: grass/trunk/imagery/i.atcorr/atmosmodel.h
===================================================================
--- grass/trunk/imagery/i.atcorr/atmosmodel.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/atmosmodel.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -39,15 +39,15 @@
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 */
Modified: grass/trunk/imagery/i.atcorr/common.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/common.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/common.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 */
Modified: grass/trunk/imagery/i.atcorr/common.h
===================================================================
--- grass/trunk/imagery/i.atcorr/common.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/common.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -38,10 +38,10 @@
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 */
@@ -48,50 +48,50 @@
/* 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 */
+ 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 */
- float usr_ph[10][83]; /* user defined model from size distribution */
- float (*ph)[10][83]; /* pointer to current active 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;
Modified: grass/trunk/imagery/i.atcorr/computations.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/computations.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/computations.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -20,22 +20,22 @@
struct OpticalAtmosProperties
{
- float rorayl, romix, roaero;
- float ddirtr, ddiftr;
- float ddirtt, ddiftt;
- float ddirta, ddifta;
- float udirtr, udiftr;
- float udirtt, udiftt;
- float udirta, udifta;
- float sphalbr, sphalbt, sphalba;
+ double rorayl, romix, roaero;
+ double ddirtr, ddiftr;
+ double ddirtt, ddiftt;
+ double ddirta, ddifta;
+ double udirtr, udiftr;
+ double udirtt, udiftt;
+ double udirta, udifta;
+ double sphalbr, sphalbt, sphalba;
};
/* To compute the molecular optical depth as a function of wavelength for any
atmosphere defined by the pressure and temperature profiles. */
-float odrayl(const AtmosModel &atms, const float wl)
+double odrayl(const AtmosModel &atms, const double wl)
{
/* air refraction index edlen 1966 / metrologia,2,71-80 putting pw=0 */
- float ak=1/wl;
+ double ak=1/wl;
double awl= wl*wl*wl*wl * 0.0254743;
double a1 = 130 - ak * ak;
double a2 = 38.9 - ak * ak;
@@ -45,12 +45,12 @@
double a = (24 * M_PI * M_PI * M_PI) * ((an * an - 1) * (an * an - 1))
* (6 + 3 * delta) / (6 - 7 * delta) / ((an * an + 2) * (an * an + 2));
- float tray = 0;
+ double tray = 0;
for(int k = 0; k < 33; k++)
{
double dppt = (288.15 / 1013.25) * (atms.p[k] / atms.t[k] + atms.p[k+1] / atms.t[k+1]) / 2;
double sr = a * dppt / awl;
- tray += (float)((atms.z[k+1] - atms.z[k]) * sr);
+ tray += (double)((atms.z[k+1] - atms.z[k]) * sr);
}
return tray;
@@ -66,13 +66,13 @@
--------
(1-w0 f)
*/
-float trunca()
+double trunca()
{
- float ptemp[83];
- float cosang[80];
- float weight[80];
- float rmu[83];
- float ga[83];
+ double ptemp[83];
+ double cosang[80];
+ double weight[80];
+ double rmu[83];
+ double ga[83];
int i;
for(i = 0; i < 83; i++) ptemp[i] = sixs_trunc.pha[i];
@@ -114,10 +114,10 @@
}
- float aa = (float)((log10(sixs_trunc.pha[kk]) - log10(sixs_trunc.pha[k])) /
+ double aa = (double)((log10(sixs_trunc.pha[kk]) - log10(sixs_trunc.pha[k])) /
(acos(rmu[kk]) - acos(rmu[k])));
- float x1 = (float)(log10(sixs_trunc.pha[kk]));
- float x2 = (float)acos(rmu[kk]);
+ double x1 = (double)(log10(sixs_trunc.pha[kk]));
+ double x2 = (double)acos(rmu[kk]);
for(i = kk + 1; i < 83; i++)
{
@@ -124,7 +124,7 @@
double a;
if(fabs(rmu[i] - 1) <= 1e-08) a = x1 - aa * x2;
else a = x1 + aa * (acos(rmu[i]) - x2);
- ptemp[i] = (float)pow(10,a);
+ ptemp[i] = (double)pow(10,a);
}
@@ -131,7 +131,7 @@
for(i = 0; i < 83; i++) sixs_trunc.pha[i] = ptemp[i];
for(i = 0; i < 80; i++) sixs_trunc.betal[i] = 0;
- float pl[83];
+ double pl[83];
#define IPL(X) ((X)+1)
@@ -138,8 +138,8 @@
for(i = 0; i < 83; i++)
{
- float x = sixs_trunc.pha[i] * ga[i];
- float rm = rmu[i];
+ double x = sixs_trunc.pha[i] * ga[i];
+ double rm = rmu[i];
pl[IPL(-1)] = 0;
pl[IPL(0)] = 1;
@@ -152,7 +152,7 @@
for(i = 0; i <= 80; i++) sixs_trunc.betal[i] *= (2 * i + 1) * 0.5f;
- float z1 = sixs_trunc.betal[0];
+ double z1 = sixs_trunc.betal[0];
for(i = 0; i <= 80; i++) sixs_trunc.betal[i] /= z1;
if(sixs_trunc.betal[80] < 0) sixs_trunc.betal[80] = 0;
@@ -168,9 +168,9 @@
unless otherwise specified by the user (using aircraft measurements).
*/
-float discre(const float ta, const float ha, const float tr, const float hr,
- const int it, const int ntd, const float yy, const float dd,
- const float ppp2, const float ppp1)
+double discre(const double ta, const double ha, const double tr, const double hr,
+ const int it, const int ntd, const double yy, const double dd,
+ const double ppp2, const double ppp1)
{
if( ha >= 7 )
{
@@ -182,14 +182,14 @@
if( it == 0 ) dt = 1e-17;
else dt = 2 * (ta + tr - yy) / (ntd - it + 1);
- float zx; /* return value */
- float ecart = 0;
+ double zx; /* return value */
+ double ecart = 0;
do {
dt = dt / 2;
double ti = yy + dt;
- float y1 = ppp2;
- float y2;
- float y3 = ppp1;
+ double y1 = ppp2;
+ double y2;
+ double y3 = ppp1;
while(true)
{
@@ -207,8 +207,8 @@
}
zx = y2;
- float cdelta = (float)(1. / (1 + ta * hr / tr / ha * exp((zx - ppp1) * (1. / hr - 1. / ha))));
- if(dd != 0) ecart = (float)fabs((dd - cdelta) / dd);
+ double cdelta = (double)(1. / (1 + ta * hr / tr / ha * exp((zx - ppp1) * (1. / hr - 1. / ha))));
+ if(dd != 0) ecart = (double)fabs((dd - cdelta) / dd);
} while((ecart > 0.75) && (it != 0));
return zx;
@@ -221,11 +221,11 @@
Compute the values of Legendre polynomials used in the successive order of
scattering method.
*/
-void kernel(const int is, float (&xpl)[2*mu + 1], float (&bp)[26][2*mu + 1], Gauss &gauss)
+void kernel(const int is, double (&xpl)[2*mu + 1], double (&bp)[26][2*mu + 1], Gauss &gauss)
{
const double rac3 = 1.7320508075688772935274463415059;
#define PSI(X) ((X)+1)
- float psl[82][2*mu + 1];
+ double psl[82][2*mu + 1];
if(is == 0)
{
@@ -238,8 +238,8 @@
double xdb = (3 * gauss.rm[STDI(j)] * gauss.rm[STDI(j)] - 1) * 0.5;
if(fabs(xdb) < 1e-30) xdb = 0;
- psl[PSI(2)][STDI(-j)] = (float)xdb;
- psl[PSI(2)][STDI(j)] = (float)xdb;
+ psl[PSI(2)][STDI(-j)] = (double)xdb;
+ psl[PSI(2)][STDI(j)] = (double)xdb;
}
psl[PSI(1)][STDI(0)] = gauss.rm[STDI(0)];
}
@@ -250,9 +250,9 @@
double x = 1 - gauss.rm[STDI(j)] * gauss.rm[STDI(j)];
psl[PSI(0)][STDI(j)] = 0;
psl[PSI(0)][STDI(-j)] = 0;
- psl[PSI(1)][STDI(-j)] = (float)sqrt(x * 0.5);
- psl[PSI(1)][STDI(j)] = (float)sqrt(x * 0.5);
- psl[PSI(2)][STDI(j)] = (float)(gauss.rm[STDI(j)] * psl[PSI(1)][STDI(j)] * rac3);
+ psl[PSI(1)][STDI(-j)] = (double)sqrt(x * 0.5);
+ psl[PSI(1)][STDI(j)] = (double)sqrt(x * 0.5);
+ psl[PSI(2)][STDI(j)] = (double)(gauss.rm[STDI(j)] * psl[PSI(1)][STDI(j)] * rac3);
psl[PSI(2)][STDI(-j)] = -psl[PSI(2)][STDI(j)];
}
@@ -270,8 +270,8 @@
psl[PSI(is - 1)][STDI(j)] = 0;
double xdb = a * pow(xx, is * 0.5);
if(fabs(xdb) < 1e-30) xdb = 0;
- psl[PSI(is)][STDI(-j)] = (float)xdb;
- psl[PSI(is)][STDI(j)] = (float)xdb;
+ psl[PSI(is)][STDI(-j)] = (double)xdb;
+ psl[PSI(is)][STDI(j)] = (double)xdb;
}
}
@@ -286,13 +286,13 @@
for(int l = k; l < ip; l++)
{
double a = (2 * l + 1.) / sqrt((l + is + 1.) * (l - is + 1.));
- double b = sqrt(float((l + is) * (l - is))) / (2. * l + 1.);
+ double b = sqrt(double((l + is) * (l - is))) / (2. * l + 1.);
for(int j = 0; j <= mu; j++)
{
double xdb = a * (gauss.rm[STDI(j)] * psl[PSI(l)][STDI(j)] - b * psl[PSI(l-1)][STDI(j)]);
if (fabs(xdb) < 1e-30) xdb = 0;
- psl[PSI(l+1)][STDI(j)] = (float)xdb;
+ psl[PSI(l+1)][STDI(j)] = (double)xdb;
if(j != 0) psl[PSI(l+1)][STDI(-j)] = ig * psl[PSI(l+1)][STDI(j)];
}
ig = -ig;
@@ -313,7 +313,7 @@
sbp += psl[PSI(l)][STDI(j)] * psl[PSI(l)][STDI(k)] * sixs_trunc.betal[l];
if(fabs(sbp) < 1e-30) sbp = 0;
- bp[j][STDI(k)] = (float)sbp;
+ bp[j][STDI(k)] = (double)sbp;
}
}
}
@@ -326,12 +326,12 @@
#define accu2 1e-3
#define mum1 (mu - 1)
-void os(const float tamoy, const float trmoy, const float pizmoy,
- const float tamoyp, const float trmoyp, float (&xl)[2*mu + 1][np],
+void os(const double tamoy, const double trmoy, const double pizmoy,
+ const double tamoyp, const double trmoyp, double (&xl)[2*mu + 1][np],
Gauss &gauss, const Altitude &alt, const GeomCond &geom)
{
- float trp = trmoy - trmoyp;
- float tap = tamoy - tamoyp;
+ double trp = trmoy - trmoyp;
+ double tap = tamoy - tamoyp;
int iplane = 0;
/* if plane observations recompute scale height for aerosol knowing:
@@ -346,16 +346,16 @@
ntp=nt-1 plane observation selected
it's a mixing rayleigh+aerosol */
- float ha = 2;
+ double ha = 2;
int snt = nt;
int ntp = snt;
if(alt.palt <= 900 && alt.palt > 0)
{
- if(tap > 1.e-03) ha = -alt.palt / (float)log(tap / tamoy);
+ if(tap > 1.e-03) ha = -alt.palt / (double)log(tap / tamoy);
ntp = snt - 1;
}
- float xmus = -gauss.rm[STDI(0)];
+ double xmus = -gauss.rm[STDI(0)];
/* compute mixing rayleigh, aerosol
case 1: pure rayleigh
@@ -362,24 +362,24 @@
case 2: pure aerosol
case 3: mixing rayleigh-aerosol */
- float h[31];
+ double h[31];
memset(h, 0, sizeof(h));
- float ch[31];
- float ydel[31];
+ double ch[31];
+ double ydel[31];
- float xdel[31];
- float altc[31];
+ double xdel[31];
+ double altc[31];
if( (tamoy <= accu2) && (trmoy > tamoy) )
{
for(int j = 0; j <= ntp; j++)
{
h[j] = j * trmoy / ntp;
- ch[j]= (float)exp(-h[j] / xmus) / 2;
+ ch[j]= (double)exp(-h[j] / xmus) / 2;
ydel[j] = 1;
xdel[j] = 0;
if (j == 0) altc[j] = 300;
- else altc[j] = -(float)log(h[j] / trmoy) * 8;
+ else altc[j] = -(double)log(h[j] / trmoy) * 8;
}
}
@@ -389,12 +389,12 @@
for(int j = 0; j <= ntp; j++)
{
h[j] = j * tamoy / ntp;
- ch[j]= (float)exp(-h[j] / xmus) / 2;
+ ch[j]= (double)exp(-h[j] / xmus) / 2;
ydel[j] = 0;
xdel[j] = pizmoy;
if (j == 0) altc[j] = 300;
- else altc[j] = -(float)log(h[j] / tamoy) * ha;
+ else altc[j] = -(double)log(h[j] / tamoy) * ha;
}
}
@@ -405,7 +405,7 @@
h[0] = 0;
ch[0] = 0.5;
altc[0] = 300;
- float zx = 300;
+ double zx = 300;
iplane = 0;
for(int it = 0; it <= ntp; it++)
@@ -414,20 +414,20 @@
else zx = discre(tamoy, ha, trmoy, 8.0, it, ntp, h[it - 1], ydel[it - 1], 300, 0);
double xx = -zx / ha;
- float ca;
+ double ca;
if( xx <= -20 ) ca = 0;
- else ca = tamoy * (float)exp(xx);
+ else ca = tamoy * (double)exp(xx);
xx = -zx / 8;
- float cr = trmoy * (float)exp(xx);
+ double cr = trmoy * (double)exp(xx);
h[it] = cr + ca;
altc[it] = zx;
- ch[it] = (float)exp(-h[it] / xmus) / 2;
+ ch[it] = (double)exp(-h[it] / xmus) / 2;
cr = cr / 8;
ca = ca / ha;
- float ratio = cr / (cr + ca);
+ double ratio = cr / (cr + ca);
xdel[it] = (1 - ratio) * pizmoy;
ydel[it] = ratio;
}
@@ -437,13 +437,13 @@
if (ntp == (snt - 1))
{
/* compute position of the plane layer */
- float taup = tap + trp;
+ double taup = tap + trp;
iplane = -1;
for(int i = 0; i <= ntp; i++) if (taup >= h[i]) iplane = i;
/* update the layer from the end to the position to update if necessary */
- float xt1 = (float)fabs(h[iplane] - taup);
- float xt2 = (float)fabs(h[iplane+1] - taup);
+ double xt1 = (double)fabs(h[iplane] - taup);
+ double xt2 = (double)fabs(h[iplane+1] - taup);
if ((xt1 > 0.0005) && (xt2 > 0.0005))
{
@@ -467,16 +467,16 @@
if ( trmoy > accu2 && tamoy > accu2 )
{
- float ca = tamoy * (float)exp(-alt.palt / ha);
- float cr = trmoy * (float)exp(-alt.palt / 8);
+ double ca = tamoy * (double)exp(-alt.palt / ha);
+ double cr = trmoy * (double)exp(-alt.palt / 8);
h[iplane] = ca + cr;
cr = cr / 8;
ca = ca / ha;
- float ratio = cr / (cr + ca);
+ double ratio = cr / (cr + ca);
xdel[iplane] = (1 - ratio) * pizmoy;
ydel[iplane] = ratio;
altc[iplane] = alt.palt;
- ch[iplane] = (float)exp(-h[iplane] / xmus) / 2;
+ ch[iplane] = (double)exp(-h[iplane] / xmus) / 2;
}
if ( trmoy > accu2 && tamoy <= accu2 )
@@ -495,28 +495,28 @@
}
- float phi = (float)geom.phirad;
+ double phi = (double)geom.phirad;
int i, k;
for(i = 0; i < np; i++) for(int m = -mu; m <= mu; m++) xl[STDI(m)][i] = 0;
/* ************ incident angle mus ******* */
- float aaaa = delta / (2 - delta);
- float ron = (1 - aaaa) / (1 + 2 * aaaa);
+ double aaaa = delta / (2 - delta);
+ double ron = (1 - aaaa) / (1 + 2 * aaaa);
/* rayleigh phase function */
- float beta0 = 1;
- float beta2 = 0.5f * ron;
+ double beta0 = 1;
+ double beta2 = 0.5f * ron;
/* fourier decomposition */
- float i1[31][2*mu + 1];
- float i2[31][2*mu + 1];
- float i3[2*mu + 1];
- float i4[2*mu + 1];
- float in[2*mu + 1];
- float inm1[2*mu + 1];
- float inm2[2*mu + 1];
+ double i1[31][2*mu + 1];
+ double i2[31][2*mu + 1];
+ double i3[2*mu + 1];
+ double i4[2*mu + 1];
+ double in[2*mu + 1];
+ double inm1[2*mu + 1];
+ double inm2[2*mu + 1];
for(i = -mu; i <= mu; i++) i4[STDI(i)] = 0;
int iborm = 80;
@@ -526,19 +526,19 @@
{
/* primary scattering */
int ig = 1;
- float roavion0 = 0;
- float roavion1 = 0;
- float roavion2 = 0;
- float roavion = 0;
+ double roavion0 = 0;
+ double roavion1 = 0;
+ double roavion2 = 0;
+ double roavion = 0;
int j;
for(j = -mu; j <= mu; j++) i3[STDI(j)] = 0;
/* kernel computations */
- float xpl[2*mu + 1];
- float bp[26][2*mu + 1];
- memset(xpl, 0, sizeof(float)*(2*mu+1));
- memset(bp, 0, sizeof(float)*26*(2*mu+1));
+ double xpl[2*mu + 1];
+ double bp[26][2*mu + 1];
+ memset(xpl, 0, sizeof(double)*(2*mu+1));
+ memset(bp, 0, sizeof(double)*26*(2*mu+1));
kernel(is,xpl,bp,gauss);
@@ -546,12 +546,12 @@
for(j = -mu; j <= mu; j++)
{
- float sa1;
- float sa2;
+ double sa1;
+ double sa2;
if((is - 2) <= 0)
{
- float spl = xpl[STDI(0)];
+ double spl = xpl[STDI(0)];
sa1 = beta0 + beta2 * xpl[STDI(j)] * spl;
sa2 = bp[0][STDI(j)];
}
@@ -564,10 +564,10 @@
for(k = 0; k <= snt; k++)
{
- float c = ch[k];
+ double c = ch[k];
double a = ydel[k];
double b = xdel[k];
- i2[k][STDI(j)] = (float)(c * (sa2 * b + sa1 * a));
+ i2[k][STDI(j)] = (double)(c * (sa2 * b + sa1 * a));
}
}
@@ -575,17 +575,17 @@
for(k = 1; k <= mu; k++)
{
i1[snt][STDI(k)] = 0;
- float zi1 = i1[snt][STDI(k)];
+ double zi1 = i1[snt][STDI(k)];
for(i = snt - 1; i >= 0; i--)
{
- float f = h[i + 1] - h[i];
+ double f = h[i + 1] - h[i];
double a = (i2[i + 1][STDI(k)] - i2[i][STDI(k)]) / f;
double b = i2[i][STDI(k)] - a * h[i];
- float c = (float)exp(-f / gauss.rm[STDI(k)]);
+ double c = (double)exp(-f / gauss.rm[STDI(k)]);
double xx = h[i] - h[i + 1] * c;
- zi1 = (float)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
+ zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
i1[i][STDI(k)] = zi1;
}
}
@@ -594,16 +594,16 @@
for(k = -mu; k <= -1; k++)
{
i1[0][STDI(k)] = 0;
- float zi1 = i1[0][STDI(k)];
+ double zi1 = i1[0][STDI(k)];
for(i = 1; i <= snt; i++)
{
- float f = h[i] - h[i - 1];
- float c = (float)exp(f / gauss.rm[STDI(k)]);
+ double f = h[i] - h[i - 1];
+ double c = (double)exp(f / gauss.rm[STDI(k)]);
double a = (i2[i][STDI(k)] -i2[i - 1][STDI(k)]) / f;
double b = i2[i][STDI(k)] - a * h[i];
double xx = h[i] - h[i - 1] * c;
- zi1 = (float)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx)/ 2);
+ zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx)/ 2);
i1[i][STDI(k)] = zi1;
}
}
@@ -659,8 +659,8 @@
if (ii2 < 1e-30) ii2 = 0;
if (ii1 < 1e-30) ii1 = 0;
- i2[i][STDI(k)] = (float)ii2;
- i2[i][STDI(-k)]= (float)ii1;
+ i2[i][STDI(k)] = (double)ii2;
+ i2[i][STDI(-k)]= (double)ii1;
}
}
}
@@ -688,8 +688,8 @@
if (ii2 < 1e-30) ii2 = 0;
if (ii1 < 1e-30) ii1 = 0;
- i2[i][STDI(k)] = (float)ii2;
- i2[i][STDI(-k)] = (float)ii1;
+ i2[i][STDI(k)] = (double)ii2;
+ i2[i][STDI(-k)] = (double)ii1;
}
}
}
@@ -699,16 +699,16 @@
for(k = 1; k <= mu; k++)
{
i1[snt][STDI(k)] = 0;
- float zi1 = i1[snt][STDI(k)];
+ double zi1 = i1[snt][STDI(k)];
for(i = snt-1; i >= 0; i--)
{
- float f = h[i + 1] - h[i];
+ double f = h[i + 1] - h[i];
double a = (i2[i + 1][STDI(k)] - i2[i][STDI(k)]) / f;
double b = i2[i][STDI(k)] - a * h[i];
- float c = (float)exp(-f / gauss.rm[STDI(k)]);
+ double c = (double)exp(-f / gauss.rm[STDI(k)]);
double xx = h[i] - h[i + 1] * c;
- zi1 = (float)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
+ zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
if (fabs(zi1) <= 1e-20) zi1 = 0;
i1[i][STDI(k)] = zi1;
}
@@ -718,16 +718,16 @@
for(k = -mu; k <= -1; k++)
{
i1[0][STDI(k)] = 0;
- float zi1 = i1[0][STDI(k)];
+ double zi1 = i1[0][STDI(k)];
for(i = 1; i <= snt; i++)
{
- float f = h[i] - h[i - 1];
- float c = (float)exp(f / gauss.rm[STDI(k)]);
+ double f = h[i] - h[i - 1];
+ double c = (double)exp(f / gauss.rm[STDI(k)]);
double a = (i2[i][STDI(k)] - i2[i - 1][STDI(k)]) / f;
double b = i2[i][STDI(k)] - a * h[i];
double xx = h[i] - h[i - 1] * c;
- zi1 = (float)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
+ zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
if (fabs(zi1) <= 1e-20) zi1 = 0;
i1[i][STDI(k)] = zi1;
@@ -745,10 +745,10 @@
/* convergence test (geometrical serie) */
if(ig > 2)
{
- float a1 = roavion2;
- float d1 = roavion1;
+ double a1 = roavion2;
+ double d1 = roavion1;
- float g1 = roavion0;
+ double g1 = roavion0;
double z = 0;
@@ -777,7 +777,7 @@
if(z < 0.0001)
{
/* successful test (geometrical serie) */
- float y1;
+ double y1;
for(int l = -mu; l <= mu; l++)
{
@@ -837,7 +837,7 @@
} while( ig <= 20 ); /* stop if order n is greater than 20 in any case */
/* sum of the fourier component s */
- float delta0s = 1;
+ double delta0s = 1;
if(is != 0) delta0s = 2;
for(k = -mu; k <= mu; k++) i4[STDI(k)] += delta0s * i3[STDI(k)];
@@ -849,8 +849,8 @@
for(int m = -mum1; m <= mum1; m++)
{
- if(m > 0) xl[STDI(m)][l] += (float)(delta0s * i3[STDI(m)] * cos(is * (phi + M_PI)));
- else xl[STDI(m)][l] += (float)(delta0s * i3[STDI(m)] * cos(is * phi));
+ if(m > 0) xl[STDI(m)][l] += (double)(delta0s * i3[STDI(m)] * cos(is * (phi + M_PI)));
+ else xl[STDI(m)][l] += (double)(delta0s * i3[STDI(m)] * cos(is * phi));
}
}
@@ -857,8 +857,8 @@
if(is == 0)
for(k = 1; k <= mum1; k++) xl[STDI(0)][0] += gauss.rm[STDI(k)] * gauss.gb[STDI(k)] * i3[STDI(-k)];
- xl[STDI(mu)][0] += (float)(delta0s * i3[STDI(mu)] * cos(is * (geom.phirad + M_PI)));
- xl[STDI(-mu)][0] += (float)(delta0s * roavion * cos(is * (geom.phirad + M_PI)));
+ xl[STDI(mu)][0] += (double)(delta0s * i3[STDI(mu)] * cos(is * (geom.phirad + M_PI)));
+ xl[STDI(-mu)][0] += (double)(delta0s * roavion * cos(is * (geom.phirad + M_PI)));
double z = 0;
for(l = -mu; l <= mu; l++)
@@ -878,8 +878,8 @@
Compute the atmospheric transmission for either a satellite or aircraft observation
as well as the spherical albedo of the atmosphere.
*/
-void iso(const float tamoy, const float trmoy, const float pizmoy,
- const float tamoyp, const float trmoyp, float (&xf)[3],
+void iso(const double tamoy, const double trmoy, const double pizmoy,
+ const double tamoyp, const double trmoyp, double (&xf)[3],
Gauss &gauss, const Altitude &alt)
{
/* molecular ratio within the layer
@@ -887,8 +887,8 @@
molecules and 2km for aerosols */
/* the optical thickness above plane are recomputed to give o.t above pla */
- float trp = trmoy - trmoyp;
- float tap = tamoy - tamoyp;
+ double trp = trmoy - trmoyp;
+ double tap = tamoy - tamoyp;
/* if plane observations recompute scale height for aerosol knowing:
the aerosol optical depth as measure from the plane = tamoyp
@@ -905,10 +905,10 @@
int snt = nt;
int iplane = 0;
int ntp = snt;
- float ha = 2.0;
+ double ha = 2.0;
if(alt.palt <= 900. && alt.palt > 0.0)
{
- if (tap > 1.e-03) ha = (float)(-alt.palt / log(tap / tamoy));
+ if (tap > 1.e-03) ha = (double)(-alt.palt / log(tap / tamoy));
else ha = 2.;
ntp = snt - 1;
}
@@ -918,10 +918,10 @@
case 2: pure aerosol
case 3: mixing rayleigh-aerosol */
- float h[31];
- float ydel[31];
- float xdel[31];
- float altc[31];
+ double h[31];
+ double ydel[31];
+ double xdel[31];
+ double altc[31];
if((tamoy <= accu2) && (trmoy > tamoy))
{
@@ -949,7 +949,7 @@
xdel[0] = 0.0;
h[0] = 0;
altc[0] = 300;
- float zx = 300;
+ double zx = 300;
iplane = 0;
for(int it = 0; it <= ntp; it++)
@@ -958,17 +958,17 @@
if (it == 0) zx = discre(tamoy,ha,trmoy,8.0,it,ntp,0,0,300.0,0.0);
else zx = discre(tamoy,ha,trmoy,8.0,it,ntp,h[it-1],ydel[it-1],300.0,0.0);
- float ca;
+ double ca;
if ((-zx / ha) < -18) ca = 0;
- else ca = (float)(tamoy * exp(-zx / ha));
+ else ca = (double)(tamoy * exp(-zx / ha));
- float cr = (float)(trmoy * exp(-zx / 8.0));
+ double cr = (double)(trmoy * exp(-zx / 8.0));
h[it] = cr + ca;
altc[it] = zx;
cr = cr / 8;
ca = ca / ha;
- float ratio = cr / (cr + ca);
+ double ratio = cr / (cr + ca);
xdel[it] = (1 - ratio) * pizmoy;
ydel[it] = ratio;
@@ -979,13 +979,13 @@
if (ntp == (snt-1))
{
/* compute position of the plane layer */
- float taup = tap + trp;
+ double taup = tap + trp;
iplane = -1;
for(int i = 0; i <= ntp; i++) if (taup >= h[i]) iplane = i;
/* update the layer from the end to the position to update if necessary */
- float xt1 = (float)fabs(h[iplane] - taup);
- float xt2 = (float)fabs(h[iplane + 1] - taup);
+ double xt1 = (double)fabs(h[iplane] - taup);
+ double xt2 = (double)fabs(h[iplane + 1] - taup);
if ((xt1 > 0.005) && (xt2 > 0.005))
{
for(int i = snt; i >= iplane + 1; i--)
@@ -1007,11 +1007,11 @@
h[iplane] = taup;
if ( trmoy > accu2 && tamoy > accu2)
{
- float ca = (float)(tamoy * exp(-alt.palt / ha));
- float cr = (float)(trmoy * exp(-alt.palt / 8.0));
+ double ca = (double)(tamoy * exp(-alt.palt / ha));
+ double cr = (double)(trmoy * exp(-alt.palt / 8.0));
cr = cr / 8;
ca = ca / ha;
- float ratio = cr / (cr + ca);
+ double ratio = cr / (cr + ca);
xdel[iplane] = (1 - ratio) * pizmoy;
ydel[iplane] = ratio;
@@ -1033,32 +1033,32 @@
}
}
- float aaaa = delta / (2-delta);
- float ron = (1 - aaaa) / (1 + 2 * aaaa);
+ double aaaa = delta / (2-delta);
+ double ron = (1 - aaaa) / (1 + 2 * aaaa);
/* rayleigh phase function */
- float beta0 = 1;
- float beta2 = 0.5f * ron;
+ double beta0 = 1;
+ double beta2 = 0.5f * ron;
/* primary scattering */
int ig = 1;
- float tavion0 = 0;
- float tavion1 = 0;
- float tavion2 = 0;
- float tavion = 0;
+ double tavion0 = 0;
+ double tavion1 = 0;
+ double tavion2 = 0;
+ double tavion = 0;
- float i1[31][2*mu + 1];
- float i2[31][2*mu + 1];
- float i3[2*mu + 1];
- float in[2*mu + 1];
- float inm1[2*mu + 1];
- float inm2[2*mu + 1];
+ double i1[31][2*mu + 1];
+ double i2[31][2*mu + 1];
+ double i3[2*mu + 1];
+ double in[2*mu + 1];
+ double inm1[2*mu + 1];
+ double inm2[2*mu + 1];
int j;
for(j = -mu; j <= mu; j++) i3[STDI(j)] = 0;
/* kernel computations */
- float xpl[2*mu + 1];
- float bp[26][2*mu + 1];
+ double xpl[2*mu + 1];
+ double bp[26][2*mu + 1];
kernel(0, xpl, bp, gauss);
for(j = -mu; j <= mu; j++)
@@ -1070,7 +1070,7 @@
{
i1[snt][STDI(k)] = 1.0;
for(int i = snt-1; i >= 0; i--)
- i1[i][STDI(k)] = (float)(exp(-(tamoy + trmoy - h[i]) / gauss.rm[STDI(k)]));
+ i1[i][STDI(k)] = (double)(exp(-(tamoy + trmoy - h[i]) / gauss.rm[STDI(k)]));
}
@@ -1111,19 +1111,19 @@
{
double ii1 = 0;
double ii2 = 0;
- float x = xdel[i];
- float y = ydel[i];
+ double x = xdel[i];
+ double y = ydel[i];
for(j = 1; j <= mu; j++)
{
- float bpjk = bp[j][STDI(k)] * x + y * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(k)]);
- float bpjmk= bp[j][STDI(-k)] * x + y * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(-k)]);
+ double bpjk = bp[j][STDI(k)] * x + y * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(k)]);
+ double bpjmk= bp[j][STDI(-k)] * x + y * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(-k)]);
ii2 += gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjk + i1[i][STDI(-j)] * bpjmk);
ii1 += gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjmk + i1[i][STDI(-j)] * bpjk);
}
- i2[i][STDI(k)] = (float)ii2;
- i2[i][STDI(-k)] = (float)ii1;
+ i2[i][STDI(k)] = (double)ii2;
+ i2[i][STDI(-k)] = (double)ii1;
}
}
@@ -1131,15 +1131,15 @@
for(k = 1; k <= mu; k++)
{
i1[snt][STDI(k)] = 0.0;
- float zi1 = i1[snt][STDI(k)];
+ double zi1 = i1[snt][STDI(k)];
for(int i = snt-1; i >= 0; i--)
{
- float f = h[i+1] - h[i];
- float a = (i2[i+1][STDI(k)] -i2[i][STDI(k)]) / f;
- float b = i2[i][STDI(k)] - a * h[i];
- float c = (float)exp(-f / gauss.rm[STDI(k)]);
- float xx = h[i] - h[i+1] * c;
+ double f = h[i+1] - h[i];
+ double a = (i2[i+1][STDI(k)] -i2[i][STDI(k)]) / f;
+ double b = i2[i][STDI(k)] - a * h[i];
+ double c = (double)exp(-f / gauss.rm[STDI(k)]);
+ double xx = h[i] - h[i+1] * c;
zi1 = c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2;
i1[i][STDI(k)] = zi1;
@@ -1150,15 +1150,15 @@
for(k = -mu; k <= -1; k++)
{
i1[0][STDI(k)] = 0;
- float zi1 = i1[0][STDI(k)];
+ double zi1 = i1[0][STDI(k)];
for(int i = 1; i <= snt; i++)
{
- float f = h[i] - h[i-1];
- float c = (float)exp(f / gauss.rm[STDI(k)]);
- float a = (i2[i][STDI(k)] - i2[i-1][STDI(k)]) / f;
- float b = i2[i][STDI(k)] - a * h[i];
- float xx = h[i] - h[i-1] * c;
+ double f = h[i] - h[i-1];
+ double c = (double)exp(f / gauss.rm[STDI(k)]);
+ double a = (i2[i][STDI(k)] - i2[i-1][STDI(k)]) / f;
+ double b = i2[i][STDI(k)] - a * h[i];
+ double xx = h[i] - h[i-1] * c;
zi1 = c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2;
i1[i][STDI(k)] = zi1;
}
@@ -1176,14 +1176,14 @@
/* convergence test (geometrical serie) */
if(ig > 2)
{
- float z = 0;
- float a1 = tavion2;
- float d1 = tavion1;
- float g1 = tavion0;
+ double z = 0;
+ double a1 = tavion2;
+ double d1 = tavion1;
+ double g1 = tavion0;
if (a1 >= accu && d1 >= accu && tavion >= accu)
{
- float y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / tavion));
- y = (float)fabs(y);
+ double y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / tavion));
+ y = (double)fabs(y);
z = z >= y ? z : y;
}
@@ -1197,8 +1197,8 @@
if(d1 == 0) continue;
if(i3[STDI(l)] == 0) continue;
- float y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / i3[STDI(l)]));
- y = (float)fabs(y);
+ double y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / i3[STDI(l)]));
+ y = (double)fabs(y);
z = z >= y ? z : y;
}
@@ -1209,7 +1209,7 @@
for(int l = -mu; l <= mu; l++)
{
if (l == 0) continue;
- float y1 = 1;
+ double y1 = 1;
d1 = inm1[STDI(l)];
g1 = in[STDI(l)];
if(d1 == 0) continue;
@@ -1224,7 +1224,7 @@
{
if (fabs(g1 - d1) >= accu)
{
- float y1 = 1 - g1 / d1;
+ double y1 = 1 - g1 / d1;
g1 = g1 / y1;
}
tavion = tavion + g1;
@@ -1248,12 +1248,12 @@
tavion = tavion + tavion0;
/* stop if order n is less than 1% of the sum */
- float z = 0;
+ double z = 0;
for(k = -mu; k <= mu; k++)
{
if(i3[STDI(k)] != 0)
{
- float y = (float)fabs(in[STDI(k)] / i3[STDI(k)]);
+ double y = (double)fabs(in[STDI(k)] / i3[STDI(k)]);
z = z >= y ? z : y;
}
}
@@ -1276,7 +1276,7 @@
To compute the atmospheric reflectance for the molecular atmosphere in
case of satellite observation.
*/
-float chand(const float xtau, const GeomCond &geom)
+double chand(const double xtau, const GeomCond &geom)
{
/* input parameters: xphi,xmus,xmuv,xtau
xphi: azimuthal difference between sun and observation (xphi=0,
@@ -1287,15 +1287,15 @@
output parameter: xrray : molecular reflectance (0.:1.)
constant : xdep: depolarization factor (0.0279) */
- const float xdep = 0.0279;
+ const double xdep = 0.0279;
- static const float as0[10] = {
+ static const double as0[10] = {
.33243832,-6.777104e-02,.16285370,1.577425e-03,-.30924818,
-1.240906e-02,-.10324388,3.241678e-02,.11493334,-3.503695e-02
};
- static const float as1[2] = { .19666292, -5.439061e-02 };
- static const float as2[2] = { .14545937,-2.910845e-02 };
+ static const double as1[2] = { .19666292, -5.439061e-02 };
+ static const double as2[2] = { .14545937,-2.910845e-02 };
double phios = (180 - geom.phi);
double xcosf1 = 1;
@@ -1344,11 +1344,11 @@
double xitot2 = xp2 + cfonc2 * fs1 * geom.xmus;
double xitot3 = xp3 + cfonc3 * fs2 * geom.xmus;
- float xrray = (float)(xitot1 * xcosf1);
- xrray += (float)(xitot2 * xcosf2 * 2);
+ double xrray = (double)(xitot1 * xcosf1);
+ xrray += (double)(xitot2 * xcosf2 * 2);
- xrray += (float)(xitot3 * xcosf3 * 2);
- xrray /= (float)geom.xmus;
+ xrray += (double)(xitot3 * xcosf3 * 2);
+ xrray /= (double)geom.xmus;
return xrray;
}
@@ -1362,13 +1362,13 @@
CHAND.f) by semi-empirical fitting of the vectorized Successive Orders of Scattering method
(Deuzé et al, 1989).
*/
-void atmref(const float tamoy, const float trmoy, const float pizmoy,
- const float tamoyp, const float trmoyp, OpticalAtmosProperties &oap,
+void atmref(const double tamoy, const double trmoy, const double pizmoy,
+ const double tamoyp, const double trmoyp, OpticalAtmosProperties &oap,
Gauss &gauss, const GeomCond &geom, const AerosolModel &aero,
const Altitude &alt)
{
- float xlm1[2 * mu + 1][np];
- float xlm2[2 * mu + 1][np];
+ double xlm1[2 * mu + 1][np];
+ double xlm2[2 * mu + 1][np];
/* atmospheric reflectances */
oap.rorayl = 0;
@@ -1377,13 +1377,13 @@
/* rayleigh reflectance 3 cases (satellite,plane,ground) */
if(alt.palt < 900 && alt.palt > 0)
{
- gauss.rm[STDI(-mu)] = -(float)geom.xmuv;
- gauss.rm[STDI(mu)] = (float)geom.xmuv;
- gauss.rm[STDI(0)] = -(float)geom.xmus;
+ gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
+ gauss.rm[STDI(mu)] = (double)geom.xmuv;
+ gauss.rm[STDI(0)] = -(double)geom.xmus;
os(0, trmoy, pizmoy, 0, trmoyp, xlm1, gauss, alt, geom);
- oap.rorayl = (float)(xlm1[STDI(-mu)][0] / geom.xmus);
+ oap.rorayl = (double)(xlm1[STDI(-mu)][0] / geom.xmus);
}
else if(alt.palt <= 0) oap.rorayl = 0;
else oap.rorayl = chand(trmoy, geom);
@@ -1399,15 +1399,15 @@
3 cases: satellite,plane,ground */
if (alt.palt > 0)
{
- gauss.rm[STDI(-mu)] = -(float)geom.xmuv;
- gauss.rm[STDI(mu)] = (float)geom.xmuv;
- gauss.rm[STDI(0)] = -(float)geom.xmus;
+ gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
+ gauss.rm[STDI(mu)] = (double)geom.xmuv;
+ gauss.rm[STDI(0)] = -(double)geom.xmus;
os(tamoy, trmoy, pizmoy, tamoyp, trmoyp, xlm2, gauss, alt, geom);
- oap.romix = (float)(xlm2[STDI(-mu)][0] / geom.xmus);
+ oap.romix = (double)(xlm2[STDI(-mu)][0] / geom.xmus);
os(tamoy, 0, pizmoy, tamoyp, 0, xlm2, gauss, alt, geom);
- oap.roaero = (float)(xlm2[STDI(-mu)][0] / geom.xmus);
+ oap.roaero = (double)(xlm2[STDI(-mu)][0] / geom.xmus);
}
else
{
@@ -1417,41 +1417,41 @@
}
-float fintexp1(const float xtau)
+double fintexp1(const double xtau)
{
/* accuracy 2e-07... for 0<xtau<1 */
- float a[6] = { -.57721566,0.99999193,-0.24991055,0.05519968,-0.00976004,0.00107857 };
- float xftau = 1;
- float xx = a[0];
+ double a[6] = { -.57721566,0.99999193,-0.24991055,0.05519968,-0.00976004,0.00107857 };
+ double xftau = 1;
+ double xx = a[0];
for(int i = 1; i <= 5; i++)
{
xftau *= xtau;
xx += a[i] * xftau;
}
- return (float)(xx - log(xtau));
+ return (double)(xx - log(xtau));
}
-float fintexp3(const float xtau)
+double fintexp3(const double xtau)
{
- return (float)((exp(-xtau) * (1. - xtau) + xtau * xtau * fintexp1(xtau)) / 2.);
+ return (double)((exp(-xtau) * (1. - xtau) + xtau * xtau * fintexp1(xtau)) / 2.);
}
/* To compute the spherical albedo of the molecular layer. */
-void csalbr(const float xtau, float& xalb)
+void csalbr(const double xtau, double& xalb)
{
- xalb = (float)((3. * xtau - fintexp3(xtau) * (4. + 2. * xtau) + 2. * exp(-xtau)));
- xalb = (float)(xalb / (4. + 3. * xtau));
+ xalb = (double)((3. * xtau - fintexp3(xtau) * (4. + 2. * xtau) + 2. * exp(-xtau)));
+ xalb = (double)(xalb / (4. + 3. * xtau));
}
-void scatra(const float taer, const float taerp,
- const float tray, const float trayp,
- const float piza, OpticalAtmosProperties& oap,
+void scatra(const double taer, const double taerp,
+ const double tray, const double trayp,
+ const double piza, OpticalAtmosProperties& oap,
Gauss &gauss, const GeomCond &geom, const Altitude &alt)
{
/* computations of the direct and diffuse transmittances
for downward and upward paths , and spherical albedo */
- float tamol,tamolp;
- float xtrans[3];
+ double tamol,tamolp;
+ double xtrans[3];
oap.ddirtt = 1; oap.ddiftt = 0;
oap.udirtt = 1; oap.udiftt = 0;
@@ -1474,13 +1474,13 @@
{
if (alt.palt > 900)
{
- oap.udiftt = (float)((2. / 3. + geom.xmuv) + (2. / 3. - geom.xmuv) * exp(-tray / geom.xmuv));
+ oap.udiftt = (double)((2. / 3. + geom.xmuv) + (2. / 3. - geom.xmuv) * exp(-tray / geom.xmuv));
- oap.udiftt = (float)(oap.udiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmuv));
- oap.ddiftt = (float)((2. / 3. + geom.xmus) + (2. / 3. - geom.xmus) * exp(-tray / geom.xmus));
- oap.ddiftt = (float)(oap.ddiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmus));
- oap.ddirtt = (float)exp(-tray / geom.xmus);
- oap.udirtt = (float)exp(-tray / geom.xmuv);
+ oap.udiftt = (double)(oap.udiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmuv));
+ oap.ddiftt = (double)((2. / 3. + geom.xmus) + (2. / 3. - geom.xmus) * exp(-tray / geom.xmus));
+ oap.ddiftt = (double)(oap.ddiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmus));
+ oap.ddirtt = (double)exp(-tray / geom.xmus);
+ oap.udirtt = (double)exp(-tray / geom.xmuv);
csalbr(tray, oap.sphalbt);
}
@@ -1488,22 +1488,22 @@
{
tamol = 0;
tamolp= 0;
- gauss.rm[STDI(-mu)] = -(float)geom.xmuv;
- gauss.rm[STDI(mu)] = (float)geom.xmuv;
- gauss.rm[STDI(0)] = (float)geom.xmus;
+ gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
+ gauss.rm[STDI(mu)] = (double)geom.xmuv;
+ gauss.rm[STDI(0)] = (double)geom.xmus;
iso(tamol, tray, piza, tamolp, trayp, xtrans, gauss, alt);
- oap.udiftt = (float)(xtrans[0] - exp(-trayp / geom.xmuv));
- oap.udirtt = (float)exp(-trayp / geom.xmuv);
- gauss.rm[STDI(-mu)] = -(float)geom.xmus;
- gauss.rm[STDI(mu)] = (float)geom.xmus;
- gauss.rm[STDI(0)] = (float)geom.xmus;
+ oap.udiftt = (double)(xtrans[0] - exp(-trayp / geom.xmuv));
+ oap.udirtt = (double)exp(-trayp / geom.xmuv);
+ gauss.rm[STDI(-mu)] = -(double)geom.xmus;
+ gauss.rm[STDI(mu)] = (double)geom.xmus;
+ gauss.rm[STDI(0)] = (double)geom.xmus;
- oap.ddiftt = (float)((2. / 3. + geom.xmus) + (2. / 3. - geom.xmus) * exp(-tray / geom.xmus));
- oap.ddiftt = (float)(oap.ddiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmus));
- oap.ddirtt = (float)exp(-tray / geom.xmus);
- oap.udirtt = (float)exp(-tray / geom.xmuv);
+ oap.ddiftt = (double)((2. / 3. + geom.xmus) + (2. / 3. - geom.xmus) * exp(-tray / geom.xmus));
+ oap.ddiftt = (double)(oap.ddiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmus));
+ oap.ddirtt = (double)exp(-tray / geom.xmus);
+ oap.udirtt = (double)exp(-tray / geom.xmuv);
csalbr(tray, oap.sphalbt);
}
@@ -1524,26 +1524,26 @@
{
tamol = 0;
tamolp= 0;
- gauss.rm[STDI(-mu)] = -(float)geom.xmuv;
- gauss.rm[STDI(mu)] = (float)geom.xmuv;
- gauss.rm[STDI(0)] = (float)geom.xmus;
+ gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
+ gauss.rm[STDI(mu)] = (double)geom.xmuv;
+ gauss.rm[STDI(0)] = (double)geom.xmus;
iso(taer, tamol, piza, taerp, tamolp, xtrans, gauss, alt);
- oap.udiftt = (float)(xtrans[0] - exp(-taerp / geom.xmuv));
- oap.udirtt = (float)exp(-taerp / geom.xmuv);
- gauss.rm[STDI(-mu)] = -(float)geom.xmus;
- gauss.rm[STDI(mu)] = (float)geom.xmus;
- gauss.rm[STDI(0)] = (float)geom.xmus;
+ oap.udiftt = (double)(xtrans[0] - exp(-taerp / geom.xmuv));
+ oap.udirtt = (double)exp(-taerp / geom.xmuv);
+ gauss.rm[STDI(-mu)] = -(double)geom.xmus;
+ gauss.rm[STDI(mu)] = (double)geom.xmus;
+ gauss.rm[STDI(0)] = (double)geom.xmus;
- float tmp_alt = alt.palt;
+ double tmp_alt = alt.palt;
alt.palt = 999;
iso(taer, tamol, piza, taerp, tamolp, xtrans, gauss, alt);
alt.palt = tmp_alt;
- oap.ddirtt = (float)exp(-taer / geom.xmus);
- oap.ddiftt = (float)(xtrans[2] - exp(-taer / geom.xmus));
- oap.sphalbt= (float)(xtrans[1] * 2);
+ oap.ddirtt = (double)exp(-taer / geom.xmus);
+ oap.ddiftt = (double)(xtrans[2] - exp(-taer / geom.xmus));
+ oap.sphalbt= (double)(xtrans[1] * 2);
if(alt.palt <= 0)
{
@@ -1559,27 +1559,27 @@
}
else if(it == 3)
{
- gauss.rm[STDI(-mu)] = -(float)geom.xmuv;
- gauss.rm[STDI(mu)] = (float)geom.xmuv;
- gauss.rm[STDI(0)] = (float)geom.xmus;
+ gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
+ gauss.rm[STDI(mu)] = (double)geom.xmuv;
+ gauss.rm[STDI(0)] = (double)geom.xmus;
iso(taer, tray, piza, taerp, trayp, xtrans, gauss, alt);
- oap.udirtt = (float)exp(-(taerp + trayp) / geom.xmuv);
- oap.udiftt = (float)(xtrans[0] - exp(-(taerp + trayp) / geom.xmuv));
- gauss.rm[STDI(-mu)] = -(float)geom.xmus;
- gauss.rm[STDI(mu)] = (float)geom.xmus;
- gauss.rm[STDI(0)] = (float)geom.xmus;
+ oap.udirtt = (double)exp(-(taerp + trayp) / geom.xmuv);
+ oap.udiftt = (double)(xtrans[0] - exp(-(taerp + trayp) / geom.xmuv));
+ gauss.rm[STDI(-mu)] = -(double)geom.xmus;
+ gauss.rm[STDI(mu)] = (double)geom.xmus;
+ gauss.rm[STDI(0)] = (double)geom.xmus;
- float tmp_alt = alt.palt;
+ double tmp_alt = alt.palt;
alt.palt = 999;
iso(taer, tray, piza, taerp, trayp, xtrans, gauss, alt);
alt.palt = tmp_alt;
- oap.ddiftt = (float)(xtrans[2] - exp(-(taer + tray) / geom.xmus));
- oap.ddirtt = (float)exp(-(taer + tray) / geom.xmus);
+ oap.ddiftt = (double)(xtrans[2] - exp(-(taer + tray) / geom.xmus));
+ oap.ddirtt = (double)exp(-(taer + tray) / geom.xmus);
- oap.sphalbt= (float)(xtrans[1] * 2);
+ oap.sphalbt= (double)(xtrans[1] * 2);
if (alt.palt <= 0)
{
@@ -1615,10 +1615,10 @@
if (((i < 9) && (sixs_disc.wldis[i] < iwave.ffu.wlinf) && (sixs_disc.wldis[i+1] < iwave.ffu.wlinf)) ||
((i > 0) && (sixs_disc.wldis[i] > iwave.ffu.wlsup) && (sixs_disc.wldis[i-1] > iwave.ffu.wlsup))) continue;
- float wl = sixs_disc.wldis[i];
+ double wl = sixs_disc.wldis[i];
/* computation of rayleigh optical depth at wl */
- float tray = odrayl(atms, wl);
- float trayp;
+ double tray = odrayl(atms, wl);
+ double trayp;
/* plane case discussed here above */
if (alt.idatmp == 0) trayp = 0;
@@ -1630,9 +1630,9 @@
/* computation of aerosol optical properties at wl */
- float taer = aerocon.taer55 * sixs_aer.ext[i] / sixs_aer.ext[3];
- float taerp = alt.taer55p * sixs_aer.ext[i] / sixs_aer.ext[3];
- float piza = sixs_aer.ome[i];
+ double taer = aerocon.taer55 * sixs_aer.ext[i] / sixs_aer.ext[3];
+ double taerp = alt.taer55p * sixs_aer.ext[i] / sixs_aer.ext[3];
+ double piza = sixs_aer.ome[i];
/* computation of atmospheric reflectances
@@ -1640,7 +1640,7 @@
roaero is aerosol ref
call plegen to decompose aerosol phase function in Betal */
- float coeff = 0;
+ double coeff = 0;
if(aero.iaer != 0)
{
for(int k = 0; k < 83; k++) sixs_trunc.pha[k] = sixs_sos.phasel[i][k];
@@ -1647,9 +1647,9 @@
coeff = trunca();
}
- float tamoy = taer * (1 - piza * coeff);
- float tamoyp = taerp * (1 - piza * coeff);
- float pizmoy = piza * (1 - coeff) / (1 - piza * coeff);
+ double tamoy = taer * (1 - piza * coeff);
+ double tamoyp = taerp * (1 - piza * coeff);
+ double pizmoy = piza * (1 - coeff) / (1 - piza * coeff);
atmref(tamoy, tray, pizmoy, tamoyp, trayp, oap, gauss, geom, aero, alt);
@@ -1685,7 +1685,7 @@
EQUIVWL.f) needed for the calculation of the downward radiation field used
in the computation of the non lambertian target contribution (main.f).
*/
-void specinterp(const float wl, float& tamoy, float& tamoyp, float& pizmoy, float& pizmoyp,
+void specinterp(const double wl, double& tamoy, double& tamoyp, double& pizmoy, double& pizmoyp,
const AerosolConcentration &aerocon, const Altitude &alt)
{
@@ -1698,28 +1698,28 @@
if(wl > sixs_disc.wldis[9]) linf = 8;
int lsup = linf + 1;
- 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];
- float alphaa = (float)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] /
+ double alphaa = (double)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] /
(sixs_aer.ext[linf] * sixs_aer.ome[linf])) / coef);
- float betaa = (float)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
- float tsca = (float)(aerocon.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));
- tamoy = (float)(aerocon.taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
- tamoyp= (float)(alt.taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+ double betaa = (double)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
+ double tsca = (double)(aerocon.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));
+ tamoy = (double)(aerocon.taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
+ tamoyp= (double)(alt.taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
pizmoy= tsca / tamoy;
pizmoyp = pizmoy;
for(int k = 0; k < 83; k++)
{
- alphaa = (float)log(sixs_sos.phasel[lsup][k] / sixs_sos.phasel[linf][k]) / coef;
- betaa = (float)(sixs_sos.phasel[linf][k] / pow(wlinf,alphaa));
- sixs_trunc.pha[k] = (float)(betaa * pow(wl,alphaa));
+ alphaa = (double)log(sixs_sos.phasel[lsup][k] / sixs_sos.phasel[linf][k]) / coef;
+ betaa = (double)(sixs_sos.phasel[linf][k] / pow(wlinf,alphaa));
+ sixs_trunc.pha[k] = (double)(betaa * pow(wl,alphaa));
}
- float coeff = trunca();
+ double coeff = trunca();
tamoy *= 1 - pizmoy * coeff;
tamoyp *= 1 - pizmoyp * coeff;
@@ -1736,36 +1736,36 @@
To compute the environment functions F(r) which allows us to account for an
inhomogeneous ground.
*/
-void enviro (const float difr, const float difa, const float r, const float palt,
- const float xmuv, float& fra, float& fae, float& fr)
+void enviro (const double difr, const double difa, const double r, const double palt,
+ const double xmuv, double& fra, double& fae, double& fr)
{
- float fae0, fra0, xlnv;
- static const float alt[16] = {
+ double fae0, fra0, xlnv;
+ static const double alt[16] = {
0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,
10.0,12.0,14.0,16.0,18.0,20.0,60.0
};
- static const float cfr1[16] = {
+ static const double cfr1[16] = {
0.730,0.710,0.656,0.606,0.560,0.516,0.473,
0.433,0.395,0.323,0.258,0.209,0.171,0.142,0.122,0.070
};
- static const float cfr2[16] = {
+ static const double cfr2[16] = {
2.8,1.51,0.845,0.634,0.524,0.465,0.429,
0.405,0.390,0.386,0.409,0.445,0.488,0.545,0.608,0.868
};
- static const float cfa1[16] = {
+ static const double cfa1[16] = {
0.239,0.396,0.588,0.626,0.612,0.505,0.454,
0.448,0.444,0.445,0.444,0.448,0.448,0.448,0.448,0.448
};
- static const float cfa2[16] = {
+ static const double cfa2[16] = {
1.40,1.20,1.02,0.86,0.74,0.56,0.46,0.42,
0.38,0.34,0.3,0.28,0.27,0.27,0.27,0.27
};
- static const float cfa3[16] = {
+ static const double cfa3[16] = {
9.17,6.26,5.48,5.16,4.74,3.65,3.24,3.15,
3.07,2.97,2.88,2.83,2.83,2.83,2.83,2.83
};
@@ -1777,26 +1777,26 @@
this calculation have been done for nadir observation
and are corrected of the effect of the view zenith angle. */
- const float a0 = 1.3347;
- const float b0 = 0.57757;
- const float a1 = -1.479;
- const float b1 = -1.5275;
+ const double a0 = 1.3347;
+ const double b0 = 0.57757;
+ const double a1 = -1.479;
+ const double b1 = -1.5275;
if (palt >= 60)
{
- fae0 = (float)(1. - 0.448 * exp( -r * 0.27) - 0.552 * exp( -r * 2.83));
- fra0 = (float)(1. - 0.930 * exp( -r * 0.080) - 0.070 * exp( -r * 1.100));
+ fae0 = (double)(1. - 0.448 * exp( -r * 0.27) - 0.552 * exp( -r * 2.83));
+ fra0 = (double)(1. - 0.930 * exp( -r * 0.080) - 0.070 * exp( -r * 1.100));
}
else
{
int i;
for(i = 0; palt >= alt[i]; i++);
- float xcfr1 = 0, xcfr2 = 0, xcfa1 = 0, xcfa2 = 0, xcfa3 = 0;
+ double xcfr1 = 0, xcfr2 = 0, xcfa1 = 0, xcfa2 = 0, xcfa3 = 0;
if ((i > 0) && (i < 16))
{
- float zmin = alt[i - 1];
- float zmax = alt[i];
+ double zmin = alt[i - 1];
+ double zmax = alt[i];
xcfr1 = cfr1[i - 1] + (cfr1[i] - cfr1[i - 1]) * (palt - zmin) / (zmax - zmin);
xcfr2 = cfr2[i - 1] + (cfr2[i] - cfr2[i - 1]) * (palt - zmin) / (zmax - zmin);
xcfa1 = cfa1[i - 1] + (cfa1[i] - cfa1[i - 1]) * (palt - zmin) / (zmax - zmin);
@@ -1813,14 +1813,14 @@
xcfa3 = cfa3[0];
}
- fra0 = (float)(1. - xcfr1 * exp(-r * xcfr2) - (1. - xcfr1) * exp(-r * 0.08));
- fae0 = (float)(1. - xcfa1 * exp(-r * xcfa2) - (1. - xcfa1) * exp(-r * xcfa3));
+ fra0 = (double)(1. - xcfr1 * exp(-r * xcfr2) - (1. - xcfr1) * exp(-r * 0.08));
+ fae0 = (double)(1. - xcfa1 * exp(-r * xcfa2) - (1. - xcfa1) * exp(-r * xcfa3));
}
/* correction of the effect of the view zenith angle */
- xlnv = (float)log(xmuv);
- fra = (float)(fra0 * (xlnv * (1 - fra0) + 1));
- fae = (float)(fae0 * ((1 + a0 * xlnv + b0 * xlnv * xlnv) + fae0 * (a1 * xlnv + b1 * xlnv * xlnv) +
+ xlnv = (double)log(xmuv);
+ fra = (double)(fra0 * (xlnv * (1 - fra0) + 1));
+ fae = (double)(fae0 * ((1 + a0 * xlnv + b0 * xlnv * xlnv) + fae0 * (a1 * xlnv + b1 * xlnv * xlnv) +
fae0 * fae0 * ((-a1 - a0) * xlnv + ( - b1 - b0) * xlnv * xlnv)));
if ((difa + difr) > 1e-03) fr = (fae * difa + fra * difr) / (difa + difr);
Modified: grass/trunk/imagery/i.atcorr/gauss.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/gauss.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/gauss.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
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,7 +42,7 @@
}
/* calculate rp & gp */
- gauss (0.f, (float)2 * M_PI, rp, gp, np);
+ gauss (0.f, (double)2 * M_PI, rp, gp, np);
}
@@ -49,7 +49,7 @@
/* 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 @@
} 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];
}
}
Modified: grass/trunk/imagery/i.atcorr/gauss.h
===================================================================
--- grass/trunk/imagery/i.atcorr/gauss.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/gauss.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -10,17 +10,17 @@
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 @@
/* 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 */
Modified: grass/trunk/imagery/i.atcorr/geomcond.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/geomcond.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/geomcond.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -46,7 +46,7 @@
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 @@
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 @@
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 @@
/* 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 @@
/* 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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
{
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 @@
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 @@
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 @@
/* 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();
}
Modified: grass/trunk/imagery/i.atcorr/geomcond.h
===================================================================
--- grass/trunk/imagery/i.atcorr/geomcond.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/geomcond.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -101,39 +101,39 @@
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();
Modified: grass/trunk/imagery/i.atcorr/interp.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/interp.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/interp.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
/* 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 @@
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 @@
}
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 @@
}
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 @@
}
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));
}
}
Modified: grass/trunk/imagery/i.atcorr/interp.h
===================================================================
--- grass/trunk/imagery/i.atcorr/interp.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/interp.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
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 */
Modified: grass/trunk/imagery/i.atcorr/iwave.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/iwave.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/iwave.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -1463,7 +1463,7 @@
}
}
-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 @@
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 @@
/* 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 @@
if (iwave > 1) {
int imax;
- float smax, sthreshold;
+ double smax, sthreshold;
imax = -1;
smax = 0;
Modified: grass/trunk/imagery/i.atcorr/iwave.h
===================================================================
--- grass/trunk/imagery/i.atcorr/iwave.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/iwave.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -42,15 +42,15 @@
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 @@
/* 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();
Modified: grass/trunk/imagery/i.atcorr/main.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/main.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/main.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -138,7 +138,7 @@
}
/* 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 @@
}
/* 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 @@
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 @@
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 @@
}
- 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 @@
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 @@
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>"),
Modified: grass/trunk/imagery/i.atcorr/transform.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/transform.cpp 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/transform.cpp 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
/* 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 @@
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 @@
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 @@
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;
*/
Modified: grass/trunk/imagery/i.atcorr/transform.h
===================================================================
--- grass/trunk/imagery/i.atcorr/transform.h 2018-08-25 17:03:59 UTC (rev 73182)
+++ grass/trunk/imagery/i.atcorr/transform.h 2018-08-25 21:42:21 UTC (rev 73183)
@@ -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 @@
/* 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 */
More information about the grass-commit
mailing list