fgmsis.cpp

6 DOF Missle Simulation

  if (flags->sw[7]) {
    double t71, t72;
    t71 = (p[11]*plg[1][2])*cd14*flags->swc[5];
    t72 = (p[12]*plg[1][2])*cd14*flags->swc[5];
    t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \
         ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \
            + t72)*stloc);
}

  /* SEMIDIURNAL */
  if (flags->sw[8]) {
    double t81, t82;
    t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
    t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
    t[7] = f2*((p[5]*plg[2][2]+ p[41]*plg[2][4] + t81)*c2tloc +(p[8]*plg[2][2] + p[42]*plg[2][4] + t82)*s2tloc);
  }

  /* TERDIURNAL */
  if (flags->sw[14]) {
    t[13] = f2 * ((p[39]*plg[3][3]+(p[93]*plg[3][4]+p[46]*plg[3][6])*cd14*flags->swc[5])* s3tloc +(p[40]*plg[3][3]+(p[94]*plg[3][4]+p[48]*plg[3][6])*cd14*flags->swc[5])* c3tloc);
}

  /* magnetic activity based on daily ap */
  if (flags->sw[9]==-1) {
    ap = input->ap_a;
    if (p[51]!=0) {
      double exp1;
      exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat))));
      if (exp1>0.99999)
        exp1=0.99999;
      if (p[24]<1.0E-4)
        p[24]=1.0E-4;
      apt[0]=sg0(exp1,p,ap->a);
      /* apt[1]=sg2(exp1,p,ap->a);
         apt[2]=sg0(exp2,p,ap->a);
         apt[3]=sg2(exp2,p,ap->a);
      */
      if (flags->sw[9]) {
        t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \
     (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \
     (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \
                 cos(hr*(tloc-p[131])));
      }
    }
  } else {
    double p44, p45;
    apd=input->ap-4.0;
    p44=p[43];
    p45=p[44];
    if (p44<0)
      p44 = 1.0E-5;
    apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44);
    if (flags->sw[9]) {
      t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \
     (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+
     (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]*
            cos(hr*(tloc-p[124])));
    }
  }

  if ((flags->sw[10])&&(input->g_long>-1000.0)) {

    /* longitudinal */
    if (flags->sw[11]) {
      t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \
     ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
      +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\
      +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \
          cos(dgtr*input->g_long) \
      +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
      +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\
      +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \
      sin(dgtr*input->g_long));
    }

    /* ut and mixed ut, longitude */
    if (flags->sw[12]){
      t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\
        (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\
        ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\
        cos(sr*(input->sec-p[71])));
      t[11]+=flags->swc[11]*\
        (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\
        cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]);
    }

    /* ut, longitude magnetic activity */
    if (flags->sw[13]) {
      if (flags->sw[9]==-1) {
        if (p[51]) {
          t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\
            ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\
             cos(dgtr*(input->g_long-p[97])))\
            +apt[0]*flags->swc[11]*flags->swc[5]*\
            (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\
            cd14*cos(dgtr*(input->g_long-p[136])) \
            +apt[0]*flags->swc[12]* \
            (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\
            cos(sr*(input->sec-p[58]));
        }
      } else {
        t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\
          ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\
          cos(dgtr*(input->g_long-p[63])))\
          +apdf*flags->swc[11]*flags->swc[5]* \
          (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \
          cd14*cos(dgtr*(input->g_long-p[118])) \
          + apdf*flags->swc[12]* \
          (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \
          cos(sr*(input->sec-p[75]));
      }
    }
  }

  /* parms not used: 82, 89, 99, 139-149 */
  tinf = p[30];
  for (i=0;i<14;i++)
    tinf = tinf + fabs(flags->sw[i+1])*t[i];
  return tinf;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

double MSIS::glob7s(double *p, struct nrlmsise_input *input,
                      struct nrlmsise_flags *flags)
{
/*    VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99
 */
  double pset=2.0;
  double t[14];
  double tt;
  double cd32, cd18, cd14, cd39;
  double p32, p18, p14, p39;
  int i,j;
  double dr=1.72142E-2;
  double dgtr=1.74533E-2;
  /* confirm parameter set */
  if (p[99]==0)
    p[99]=pset;
  if (p[99]!=pset) {
    printf("Wrong parameter set for glob7s\n");
    return -1;
  }
  for (j=0;j<14;j++)
    t[j]=0.0;
  cd32 = cos(dr*(input->doy-p[31]));
  cd18 = cos(2.0*dr*(input->doy-p[17]));
  cd14 = cos(dr*(input->doy-p[13]));
  cd39 = cos(2.0*dr*(input->doy-p[38]));
  p32=p[31];
  p18=p[17];
  p14=p[13];
  p39=p[38];

  /* F10.7 */
  t[0] = p[21]*dfa;

  /* time independent */
  t[1]=p[1]*plg[0][2] + p[2]*plg[0][4] + p[22]*plg[0][6] + p[26]*plg[0][1] + p[14]*plg[0][3] + p[59]*plg[0][5];

        /* SYMMETRICAL ANNUAL */
  t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32;

        /* SYMMETRICAL SEMIANNUAL */
  t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18;

        /* ASYMMETRICAL ANNUAL */
  t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14;

  /* ASYMMETRICAL SEMIANNUAL */
  t[5]=(p[37]*plg[0][1])*cd39;

        /* DIURNAL */
  if (flags->sw[7]) {
    double t71, t72;
    t71 = p[11]*plg[1][2]*cd14*flags->swc[5];
    t72 = p[12]*plg[1][2]*cd14*flags->swc[5];
    t[6] = ((p[3]*plg[1][1] + p[4]*plg[1][3] + t71) * ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + t72) * stloc) ;
  }

  /* SEMIDIURNAL */
  if (flags->sw[8]) {
    double t81, t82;
    t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
    t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
    t[7] = ((p[5]*plg[2][2] + p[41]*plg[2][4] + t81) * c2tloc + (p[8]*plg[2][2] + p[42]*plg[2][4] + t82) * s2tloc);
  }

  /* TERDIURNAL */
  if (flags->sw[14]) {
    t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
  }

  /* MAGNETIC ACTIVITY */
  if (flags->sw[9]) {
    if (flags->sw[9]==1)
      t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]);
    if (flags->sw[9]==-1)
      t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]);
  }

  /* LONGITUDINAL */
  if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0))) {
    t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*cos(dr*(input->doy-p[81]))\
            +p[85]*flags->swc[6]*cos(2.0*dr*(input->doy-p[86])))\
      +p[83]*flags->swc[3]*cos(dr*(input->doy-p[84]))\
      +p[87]*flags->swc[4]*cos(2.0*dr*(input->doy-p[88])))\
      *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
      +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\
      )*cos(dgtr*input->g_long)\
      +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
      +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\
      )*sin(dgtr*input->g_long));
  }
  tt=0;
  for (i=0;i<14;i++)
    tt+=fabs(flags->sw[i+1])*t[i];
  return tt;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void MSIS::gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
                  struct nrlmsise_output *output)
{
  double xlat;
  double xmm;
  int mn3 = 5;
  double zn3[5]={32.5,20.0,15.0,10.0,0.0};
  int mn2 = 4;
  double zn2[4]={72.5,55.0,45.0,32.5};
  double altt;
  double zmix=62.5;
  double tmp;
  double dm28m;
  double tz;
  double dmc;
  double dmr;
  double dz28;
  struct nrlmsise_output soutput;
  int i;

  tselec(flags);

  /* Latitude variation of gravity (none for sw[2]=0) */
  xlat=input->g_lat;
  if (flags->sw[2]==0)
    xlat=45.0;
  glatf(xlat, &gsurf, &re);

  xmm = pdm[2][4];

  /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
  if (input->alt>zn2[0])
    altt=input->alt;
  else
    altt=zn2[0];

  tmp=input->alt;
  input->alt=altt;
  gts7(input, flags, &soutput);
  altt=input->alt;
  input->alt=tmp;
  if (flags->sw[0])   /* metric adjustment */
    dm28m=dm28*1.0E6;
  else
    dm28m=dm28;
  output->t[0]=soutput.t[0];
  output->t[1]=soutput.t[1];
  if (input->alt>=zn2[0]) {
    for (i=0;i<9;i++)
      output->d[i]=soutput.d[i];
    return;
  }

/*       LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
 *         Temperature at nodes and gradients at end nodes
 *         Inverse temperature a linear function of spherical harmonics
 */
  meso_tgn2[0]=meso_tgn1[1];
  meso_tn2[0]=meso_tn1[4];
        meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags));
        meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags));
        meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags));
  meso_tgn2[1]=pavgm[8]*pma[9][0]*(1.0+flags->sw[20]*flags->sw[22]*glob7s(pma[9], input, flags))*meso_tn2[3]*meso_tn2[3]/(pow((pma[2][0]*pavgm[2]),2.0));
  meso_tn3[0]=meso_tn2[3];

  if (input->alt<zn3[0]) {
/*       LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
 *         Temperature at nodes and gradients at end nodes
 *         Inverse temperature a linear function of spherical harmonics
 */
    meso_tgn3[0]=meso_tgn2[1];
    meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags));
    meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags));
    meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags));
    meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags));
    meso_tgn3[1]=pma[7][0]*pavgm[7]*(1.0+flags->sw[22]*glob7s(pma[7], input, flags)) *meso_tn3[4]*meso_tn3[4]/(pow((pma[6][0]*pavgm[6]),2.0));
  }

        /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */

  dmc=0;
  if (input->alt>zmix)
    dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
  dz28=soutput.d[2];

  /**** N2 density ****/
  dmr=soutput.d[2] / dm28m - 1.0;
  output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
  output->d[2]=output->d[2] * (1.0 + dmr*dmc);

  /**** HE density ****/
  dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
  output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc);

  /**** O density ****/
  output->d[1] = 0;
  output->d[8] = 0;

  /**** O2 density ****/
  dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
  output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc);

  /**** AR density ***/
  dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
  output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc);

  /**** Hydrogen density ****/
  output->d[6] = 0;

  /**** Atomic nitrogen density ****/
  output->d[7] = 0;

  /**** Total mass density */
  output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] +
                     28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4]
                     + output->d[6] + 14.0 * output->d[7]);

  if (flags->sw[0])
    output->d[5]=output->d[5]/1000;

  /**** temperature at altitude ****/
  dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3,
                   mn2, zn2, meso_tn2, meso_tgn2);
  output->t[1]=tz;

}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void MSIS::gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
                   struct nrlmsise_output *output)
{
  gtd7(input, flags, output);
  output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] +
                   28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4]
                   + output->d[6] + 14.0 * output->d[7] + 16.0 * output->d[8]);
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void MSIS::ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
                  struct nrlmsise_output *output, double press)
{
  double bm = 1.3806E-19;
  double rgas = 831.4;
  double test = 0.00043;
  double ltest = 12;
  double pl, p;
  double zi = 0.0;
  double z;
  double cl, cl2;
  double ca, cd;
  double xn, xm, diff;
  double g, sh;
  int l;
  pl = log10(press);
  if (pl >= -5.0) {
    if (pl>2.5)
      zi = 18.06 * (3.00 - pl);
    else if ((pl>0.075) && (pl<=2.5))
      zi = 14.98 * (3.08 - pl);
    else if ((pl>-1) && (pl<=0.075))
      zi = 17.80 * (2.72 - pl);
    else if ((pl>-2) && (pl<=-1))
      zi = 14.28 * (3.64 - pl);
    else if ((pl>-4) && (pl<=-2))
      zi = 12.72 * (4.32 -pl);
    else if (pl<=-4)
      zi = 25.3 * (0.11 - pl);
    cl = input->g_lat/90.0;
    cl2 = cl*cl;
    if (input->doy<182)
      cd = (1.0 - (double) input->doy) / 91.25;
    else
      cd = ((double) input->doy) / 91.25 - 3.0;
    ca = 0;
    if ((pl > -1.11) && (pl<=-0.23))
      ca = 1.0;
    if (pl > -0.23)
      ca = (2.79 - pl) / (2.79 + 0.23);
    if ((pl <= -1.11) && (pl>-3))
      ca = (-2.93 - pl)/(-2.93 + 1.11);
    z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl;
  } else
    z = 22.0 * pow((pl + 4.0),2.0) + 110.0;

  /* iteration  loop */
  l = 0;
  do {
    l++;
    input->alt = z;
    gtd7(input, flags, output);
    z = input->alt;
    xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7];
    p = bm * xn * output->t[1];