mirror of
https://github.com/eddyem/BTA_utils.git
synced 2025-12-06 10:45:14 +03:00
76 lines
3.2 KiB
C
76 lines
3.2 KiB
C
static int main(int argc, char *argv[]) __attribute((unused));
|
|
static void usage(void) __attribute((unused));
|
|
#define inline
|
|
#include "airmass.inc"
|
|
#undef inline
|
|
|
|
#include "bta_site.h"
|
|
|
|
void calc_airmass(
|
|
// in
|
|
double daynum, // Day number from beginning of year, 0 = midnight Jan 1
|
|
double h0,// Relative humidity in percent
|
|
double p0, // local pressure in mmHg
|
|
double t0, // temperature in kelvins
|
|
double z, // zenith distance in degr
|
|
// out
|
|
double *am, // AIRMASS
|
|
double *acd, // column density
|
|
double *wcd, // water vapor column density
|
|
double *wam // water vapor airmass
|
|
){
|
|
int i;
|
|
double colz, col0, cosyearfrac, sinz, rmsl, bigrmsl,
|
|
rjunk, coslatsq, lambda=5500., alt=TELALT,
|
|
lat=TELLAT, betapoly[5]={-0.0065107, -4.5403e-06,
|
|
3.6599e-07, -2.2174e-09, 7.9392e-12}, r1poly[11]={17.204, 8.9155e-03,
|
|
-3.6420e-03, 2.5617e-05, 2.4796e-07, -1.2774e-08, 1.3017e-10,
|
|
2.0151e-12, -2.6985e-14, -1.0397e-16, 1.4849e-18}, hvec[(NTLEVELS+1)]=
|
|
{0., 11000., 20000., 32000., 47000., 52000., 61000., 79000., 88743.};
|
|
p0 *= 101325. / 760.;
|
|
relhumid = h0 / 100.;
|
|
sinz=sin(z*DEG2RAD);
|
|
cee=CONSTC*(2.87566E-4+134.12/(lambda*lambda)+3.777E8*pow(lambda, -4.));
|
|
cosyearfrac=cos((daynum-202.)*(360.*DEG2RAD/365.));
|
|
coslatsq=cos(lat*DEG2RAD);
|
|
coslatsq*=coslatsq;
|
|
rmsl=sqrt(((POLAR*POLAR*POLAR*POLAR)+(EQUATOR*EQUATOR*EQUATOR*EQUATOR-
|
|
POLAR*POLAR*POLAR*POLAR)*coslatsq)/((POLAR*POLAR)+(EQUATOR*EQUATOR-
|
|
POLAR*POLAR)*coslatsq));
|
|
bigrmsl=(-STANDARDR02)/rmsl;
|
|
/* Calculate bigrvec, the bigr at the bottom of each temperature layer */
|
|
*bigrvec=(-STANDARDR02)/(rmsl+alt);
|
|
rjunk=r1poly[10];
|
|
for (i=9;i>=0;i--) rjunk=rjunk*lat+((i%2) ? r1poly[i]*cosyearfrac :
|
|
r1poly[i]);
|
|
bigrvec[1]=(-STANDARDR02)/(rmsl+rjunk*1000.);
|
|
for (i=2;i<(NTLEVELS+1);i++) bigrvec[i]=hvec[i]+bigrmsl;
|
|
|
|
/* Set up our temperature profile for the troposphere/lower stratosphere */
|
|
*betavec=betapoly[4];
|
|
for (i=3;i>=0;i--) *betavec=*betavec*lat+((i%2) ? betapoly[i]*cosyearfrac :
|
|
betapoly[i]);
|
|
*betavec*=(-rmsl/bigrvec[1]);
|
|
*tvec=t0;
|
|
tvec[1]=t0+*betavec*(bigrvec[1]-*bigrvec);
|
|
betavec[1]=(tvec[2]-tvec[1])/(bigrvec[2]-bigrvec[1]);
|
|
|
|
/* Compute the density at the bottom of each temperature layer */
|
|
*rhovec=p0/(AIRCONSTANT*t0);
|
|
for (i=0;i<(NTLEVELS-1);i++) rhovec[i+1]=frho(bigrvec[i+1]-bigrvec[i],
|
|
rhovec[i], tvec[i], betavec[i]);
|
|
|
|
const6=musquare(*rhovec)*sinz*sinz/(*bigrvec*(*bigrvec)); /* for z */
|
|
colz=qromb(columndensityint, *bigrvec, bigrvec[NTLEVELS], 1.E-8);
|
|
const6=0.; /* equivalent to setting z = 0. */
|
|
col0=qromb(columndensityint, *bigrvec, bigrvec[NTLEVELS], 1.E-8);
|
|
*am = colz/col0; *acd = 0.1*colz;
|
|
if(h0 < 0.1){ *wcd = 0.; *wam = 0.; return;}
|
|
const7=relhumid/(bigrvec[2]-bigrvec[1]);
|
|
const6=musquare(*rhovec)*sinz*sinz/(*bigrvec*(*bigrvec)); /* for z */
|
|
colz=qromb(vaporcolumndensityint, *bigrvec, bigrvec[2], 1.E-8);
|
|
const6=0.; /* equivalent to setting z = 0. */
|
|
col0=qromb(vaporcolumndensityint, *bigrvec, bigrvec[2], 1.E-8);
|
|
*wcd = 0.1*colz; *wam = colz/col0;
|
|
}
|