view release on metacpan or  search on metacpan

palsrc/palRefro.c  view on Meta::CPAN

  psat = pow(10.0,(0.7859+0.03477*tdc)/(1.0+0.00412*tdc)) *
    (1.0+pmbok*(4.5e-6+6.0e-10*tdc*tdc));
  if (pmbok > 0.0) {
    pwo = rhok*psat/(1.0-(1.0-rhok)*psat/pmbok);
  } else {
    pwo = 0.0;
  }
  w = pwo*(1.0-DMW/DMD)*gamma/(DELTA-gamma);
  c1 = a*(pmbok+w)/tdkok;
  if (optic) {
    c2 = (a*w+11.2684e-6*pwo)/tdkok;
  } else {
    c2 = (a*w+6.3938e-6*pwo)/tdkok;
  }
  c3 = (gamma-1.0)*alpha*c1/tdkok;
  c4 = (DELTA-1.0)*alpha*c2/tdkok;
  if (optic) {
    c5 = 0.0;
    c6 = 0.0;
  } else {
    c5 = 375463e-6*pwo/tdkok;
    c6 = c5*delm2*alpha/(tdkok*tdkok);
  }

  /*  conditions at the observer. */
  r0 = S+hmok;
  pal1Atmt(r0,tdkok,alpha,gamm2,delm2,c1,c2,c3,c4,c5,c6,
           r0,&tempo,&dn0,&rdndr0);
  sk0 = dn0*r0*sin(zobs2);
  f0 = refi(dn0,rdndr0);

  /*  conditions in the troposphere at the tropopause. */
  rt = S+DMAX(HT,hmok);
  pal1Atmt(r0,tdkok,alpha,gamm2,delm2,c1,c2,c3,c4,c5,c6,
           rt,&tt,&dnt,&rdndrt);
  sine = sk0/(rt*dnt);
  zt = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
  ft = refi(dnt,rdndrt);

  /*  conditions in the stratosphere at the tropopause. */
  pal1Atms(rt,tt,dnt,gamal,rt,&dnts,&rdndrp);
  sine = sk0/(rt*dnts);
  zts = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
  fts = refi(dnts,rdndrp);

  /*  conditions at the stratosphere limit. */
  rs = S+HS;
  pal1Atms(rt,tt,dnt,gamal,rs,&dns,&rdndrs);
  sine = sk0/(rs*dns);
  zs = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
  fs = refi(dns,rdndrs);

  /*  variable initialization to avoid compiler warning. */
  reft = 0.0;

  /*  integrate the refraction integral in two parts;  first in the
   *  troposphere (k=1), then in the stratosphere (k=2). */

  for (k=1; k<=2; k++) {

    /*     initialize previous refraction to ensure at least two iterations. */
    refold = 1.0;

    /*     start off with 8 strips. */
    is = 8;

    /*     start z, z range, and start and end values. */
    if (k==1) {
      z0 = zobs2;
      zrange = zt-z0;
      fb = f0;
      ff = ft;
    } else {
      z0 = zts;
      zrange = zs-z0;
      fb = fts;
      ff = fs;
    }

    /*     sums of odd and even values. */
    fo = 0.0;
    fe = 0.0;

    /*     first time through the loop we have to do every point. */
    n = 1;

    /*     start of iteration loop (terminates at specified precision). */
    loop = 1;
    while (loop) {

      /*        strip width. */
      h = zrange/((double)is);

      /*        initialize distance from earth centre for quadrature pass. */
      if (k == 1) {
        r = r0;
      } else {
        r = rt;
      }

      /*        one pass (no need to compute evens after first time). */
      for (i=1; i<is; i+=n) {

              /*           sine of observed zenith distance. */
        sz = sin(z0+h*(double)(i));

        /*           find r (to the nearest metre, maximum four iterations). */
        if (sz > 1e-20) {
          w = sk0/sz;
          rg = r;
          dr = 1.0e6;
          j = 0;
          while ( fabs(dr) > 1.0 && j < 4 ) {
            j++;
            if (k==1) {
              pal1Atmt(r0,tdkok,alpha,gamm2,delm2,
                       c1,c2,c3,c4,c5,c6,rg,&tg,&dn,&rdndr);
            } else {
              pal1Atms(rt,tt,dnt,gamal,rg,&dn,&rdndr);
            }
            dr = (rg*dn-w)/(dn+rdndr);
            rg = rg-dr;
          }
          r = rg;
        }

        /*           find the refractive index and integrand at r. */
        if (k==1) {
          pal1Atmt(r0,tdkok,alpha,gamm2,delm2,
                   c1,c2,c3,c4,c5,c6,r,&t,&dn,&rdndr);
        } else {
          pal1Atms(rt,tt,dnt,gamal,r,&dn,&rdndr);
        }
        f = refi(dn,rdndr);

        /*           accumulate odd and (first time only) even values. */
        if (n==1 && i%2 == 0) {
          fe += f;
        } else {
          fo += f;
        }
      }

      /*        evaluate the integrand using simpson's rule. */
      refp = h*(fb+4.0*fo+2.0*fe+ff)/3.0;

      /*        has the required precision been achieved (or can't be)? */
      if (fabs(refp-refold) > tol && is < ISMAX) {

        /*           no: prepare for next iteration.*/

        /*           save current value for convergence test. */
        refold = refp;

        /*           double the number of strips. */
        is += is;

        /*           sum of all current values = sum of next pass's even values. */
        fe += fo;

        /*           prepare for new odd values. */
        fo = 0.0;

        /*           skip even values next time. */
        n = 2;
      } else {