fgmsis.cpp

6 DOF Missle Simulation

    if (flags->sw[0])
      p = p*1.0E-6;
    diff = pl - log10(p);
    if (sqrt(diff*diff)<test)
      return;
    if (l==ltest) {
      printf("ERROR: ghp7 not converging for press %e, diff %e",press,diff);
      return;
    }
    xm = output->d[5] / xn / 1.66E-24;
    if (flags->sw[0])
      xm = xm * 1.0E3;
    g = gsurf / (pow((1.0 + z/re),2.0));
    sh = rgas * output->t[1] / (xm * g);

    /* new altitude estimate using scale height */
    if (l <  6)
      z = z - sh * diff * 2.302;
    else
      z = z - sh * diff;
  } while (1==1);
}

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

void MSIS::gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
                  struct nrlmsise_output *output)
{
/*     Thermospheric portion of NRLMSISE-00
 *     See GTD7 for more extensive comments
 *     alt > 72.5 km!
 */
  double za;
  int i, j;
  double ddum, z;
  double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
  double tinf;
  int mn1 = 5;
  double g0;
  double tlb;
  double s, z0, t0, tr12;
  double db01, db04, db14, db16, db28, db32, db40, db48;
  double zh28, zh04, zh16, zh32, zh40, zh01, zh14;
  double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14;
  double xmd;
  double b28, b04, b16, b32, b40, b01, b14;
  double tz;
  double g28, g4, g16, g32, g40, g1, g14;
  double zhf, xmm;
  double zc04, zc16, zc32, zc40, zc01, zc14;
  double hc04, hc16, hc32, hc40, hc01, hc14;
  double hcc16, hcc32, hcc01, hcc14;
  double zcc16, zcc32, zcc01, zcc14;
  double rc16, rc32, rc01, rc14;
  double rl;
  double g16h, db16h, tho, zsht, zmho, zsho;
  double dgtr=1.74533E-2;
  double dr=1.72142E-2;
  double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
  double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
  double dd;
  double hc216, hcc232;
  za = pdl[1][15];
  zn1[0] = za;
  for (j=0;j<9;j++)
    output->d[j]=0;

  /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
  if (input->alt>zn1[0])
    tinf = ptm[0]*pt[0] * \
      (1.0+flags->sw[16]*globe7(pt,input,flags));
  else
    tinf = ptm[0]*pt[0];
  output->t[0]=tinf;

  /*  GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
  if (input->alt>zn1[4])
    g0 = ptm[3]*ps[0] * \
      (1.0+flags->sw[19]*globe7(ps,input,flags));
  else
    g0 = ptm[3]*ps[0];
  tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
  s = g0 / (tinf - tlb);

/*      Lower thermosphere temp variations not significant for
 *       density above 300 km */
  if (input->alt<300.0) {
    meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags));
    meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags));
    meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags));
    meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags));
    meso_tgn1[1]=ptm[8]*pma[8][0]*(1.0+flags->sw[18]*flags->sw[20]*glob7s(pma[8], input, flags))*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
  } else {
    meso_tn1[1]=ptm[6]*ptl[0][0];
    meso_tn1[2]=ptm[2]*ptl[1][0];
    meso_tn1[3]=ptm[7]*ptl[2][0];
    meso_tn1[4]=ptm[4]*ptl[3][0];
    meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
  }

  z0 = zn1[3];
  t0 = meso_tn1[3];
  tr12 = 1.0;

  /* N2 variation factor at Zlb */
  g28=flags->sw[21]*globe7(pd[2], input, flags);

  /* VARIATION OF TURBOPAUSE HEIGHT */
  zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13])));
  output->t[0]=tinf;
  xmm = pdm[2][4];
  z = input->alt;


        /**** N2 DENSITY ****/

  /* Diffusive density at Zlb */
  db28 = pdm[2][0]*exp(g28)*pd[2][0];
  /* Diffusive density at Alt */
  output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
  dd=output->d[2];
  /* Turbopause */
  zh28=pdm[2][2]*zhf;
  zhm28=pdm[2][3]*pdl[1][5];
  xmd=28.0-xmm;
  /* Mixed density at Zlb */
  b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
  if ((flags->sw[15])&&(z<=altl[2])) {
    /*  Mixed density at Alt */
    dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    /*  Net density at Alt */
    output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
  }


        /**** HE DENSITY ****/

  /*   Density variation factor at Zlb */
  g4 = flags->sw[21]*globe7(pd[0], input, flags);
  /*  Diffusive density at Zlb */
  db04 = pdm[0][0]*exp(g4)*pd[0][0];
        /*  Diffusive density at Alt */
  output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
  dd=output->d[0];
  if ((flags->sw[15]) && (z<altl[0])) {
    /*  Turbopause */
    zh04=pdm[0][2];
    /*  Mixed density at Zlb */
    b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    /*  Mixed density at Alt */
    dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    zhm04=zhm28;
    /*  Net density at Alt */
    output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
    /*  Correction to specified mixing ratio at ground */
    rl=log(b28*pdm[0][1]/b04);
    zc04=pdm[0][4]*pdl[1][0];
    hc04=pdm[0][5]*pdl[1][1];
    /*  Net density corrected at Alt */
    output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
  }


        /**** O DENSITY ****/

  /*  Density variation factor at Zlb */
  g16= flags->sw[21]*globe7(pd[1],input,flags);
  /*  Diffusive density at Zlb */
  db16 =  pdm[1][0]*exp(g16)*pd[1][0];
        /*   Diffusive density at Alt */
  output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
  dd=output->d[1];
  if ((flags->sw[15]) && (z<=altl[1])) {
    /*   Turbopause */
    zh16=pdm[1][2];
    /*  Mixed density at Zlb */
    b16=densu(zh16,db16,tinf,tlb,16.0-xmm,(alpha[1]-1.0), &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    /*  Mixed density at Alt */
    dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    zhm16=zhm28;
    /*  Net density at Alt */
    output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.);
    rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0));
    hc16=pdm[1][5]*pdl[1][3];
    zc16=pdm[1][4]*pdl[1][2];
    hc216=pdm[1][5]*pdl[1][4];
    output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216);
    /*   Chemistry correction */
    hcc16=pdm[1][7]*pdl[1][13];
    zcc16=pdm[1][6]*pdl[1][12];
    rc16=pdm[1][3]*pdl[1][14];
    /*  Net density corrected at Alt */
    output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
  }


        /**** O2 DENSITY ****/

        /*   Density variation factor at Zlb */
  g32= flags->sw[21]*globe7(pd[4], input, flags);
        /*  Diffusive density at Zlb */
  db32 = pdm[3][0]*exp(g32)*pd[4][0];
        /*   Diffusive density at Alt */
  output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
  dd=output->d[3];
  if (flags->sw[15]) {
    if (z<=altl[3]) {
      /*   Turbopause */
      zh32=pdm[3][2];
      /*  Mixed density at Zlb */
      b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      /*  Mixed density at Alt */
      dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
      zhm32=zhm28;
      /*  Net density at Alt */
      output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
      /*   Correction to specified mixing ratio at ground */
      rl=log(b28*pdm[3][1]/b32);
      hc32=pdm[3][5]*pdl[1][7];
      zc32=pdm[3][4]*pdl[1][6];
      output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32);
    }
    /*  Correction for general departure from diffusive equilibrium above Zlb */
    hcc32=pdm[3][7]*pdl[1][22];
    hcc232=pdm[3][7]*pdl[0][22];
    zcc32=pdm[3][6]*pdl[1][21];
    rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.));
    /*  Net density corrected at Alt */
    output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
  }


        /**** AR DENSITY ****/

        /*   Density variation factor at Zlb */
  g40= flags->sw[20]*globe7(pd[5],input,flags);
        /*  Diffusive density at Zlb */
  db40 = pdm[4][0]*exp(g40)*pd[5][0];
  /*   Diffusive density at Alt */
  output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
  dd=output->d[4];
  if ((flags->sw[15]) && (z<=altl[4])) {
    /*   Turbopause */
    zh40=pdm[4][2];
    /*  Mixed density at Zlb */
    b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    /*  Mixed density at Alt */
    dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    zhm40=zhm28;
    /*  Net density at Alt */
    output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
    /*   Correction to specified mixing ratio at ground */
    rl=log(b28*pdm[4][1]/b40);
    hc40=pdm[4][5]*pdl[1][9];
    zc40=pdm[4][4]*pdl[1][8];
    /*  Net density corrected at Alt */
    output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
    }


        /**** HYDROGEN DENSITY ****/

        /*   Density variation factor at Zlb */
  g1 = flags->sw[21]*globe7(pd[6], input, flags);
        /*  Diffusive density at Zlb */
  db01 = pdm[5][0]*exp(g1)*pd[6][0];
        /*   Diffusive density at Alt */
  output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
  dd=output->d[6];
  if ((flags->sw[15]) && (z<=altl[6])) {
    /*   Turbopause */
    zh01=pdm[5][2];
    /*  Mixed density at Zlb */
    b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    /*  Mixed density at Alt */
    dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    zhm01=zhm28;
    /*  Net density at Alt */
    output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
    /*   Correction to specified mixing ratio at ground */
    rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01);
    hc01=pdm[5][5]*pdl[1][11];
    zc01=pdm[5][4]*pdl[1][10];
    output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01);
    /*   Chemistry correction */
    hcc01=pdm[5][7]*pdl[1][19];
    zcc01=pdm[5][6]*pdl[1][18];
    rc01=pdm[5][3]*pdl[1][20];
    /*  Net density corrected at Alt */
    output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
}


        /**** ATOMIC NITROGEN DENSITY ****/

  /*   Density variation factor at Zlb */
  g14 = flags->sw[21]*globe7(pd[7],input,flags);
        /*  Diffusive density at Zlb */
  db14 = pdm[6][0]*exp(g14)*pd[7][0];
        /*   Diffusive density at Alt */
  output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
  dd=output->d[7];
  if ((flags->sw[15]) && (z<=altl[7])) {
    /*   Turbopause */
    zh14=pdm[6][2];
    /*  Mixed density at Zlb */
    b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    /*  Mixed density at Alt */
    dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
    zhm14=zhm28;
    /*  Net density at Alt */
    output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
    /*   Correction to specified mixing ratio at ground */
    rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14);
    hc14=pdm[6][5]*pdl[0][1];
    zc14=pdm[6][4]*pdl[0][0];
    output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14);
    /*   Chemistry correction */
    hcc14=pdm[6][7]*pdl[0][4];
    zcc14=pdm[6][6]*pdl[0][3];
    rc14=pdm[6][3]*pdl[0][5];
    /*  Net density corrected at Alt */
    output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
  }


        /**** Anomalous OXYGEN DENSITY ****/

  g16h = flags->sw[21]*globe7(pd[8],input,flags);
  db16h = pdm[7][0]*exp(g16h)*pd[8][0];
  tho = pdm[7][9]*pdl[0][6];
  dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
  zsht=pdm[7][5];
  zmho=pdm[7][4];
  zsho=scalh(zmho,16.0,tho);
  output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));


  /* 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]);
  db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);



  /* temperature */
  z = sqrt(input->alt*input->alt);
  ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
  if (flags->sw[0]) {
    for(i=0;i<9;i++)
      output->d[i]=output->d[i]*1.0E6;
    output->d[5]=output->d[5]/1000;
  }
}


//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//    The bitmasked value choices are as follows:
//    unset: In this case (the default) JSBSim would only print
//       out the normally expected messages, essentially echoing
//       the config files as they are read. If the environment
//       variable is not set, debug_lvl is set to 1 internally
//    0: This requests JSBSim not to output any messages
//       whatsoever.
//    1: This value explicity requests the normal JSBSim
//       startup messages
//    2: This value asks for a message to be printed out when
//       a class is instantiated
//    4: When this value is set, a message is displayed when a
//       FGModel object executes its Run() method
//    8: When this value is set, various runtime state variables
//       are printed out periodically
//    16: When set various parameters are sanity checked and
//       a message is printed out when they go out of bounds

void MSIS::Debug(int from)
{
  if (debug_lvl <= 0) return;

  if (debug_lvl & 1) { // Standard console startup message output
    if (from == 0) { // Constructor
    }
  }
  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
    if (from == 0) cout << "Instantiated: MSIS" << endl;
    if (from == 1) cout << "Destroyed:    MSIS" << endl;
  }
  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
  }
  if (debug_lvl & 8 ) { // Runtime state variables
  }
  if (debug_lvl & 16) { // Sanity checking
  }
  if (debug_lvl & 32) { // Turbulence
    if (first_pass && from == 2) {
      cout << "vTurbulenceNED(X), vTurbulenceNED(Y), vTurbulenceNED(Z), "
           << "vTurbulenceGrad(X), vTurbulenceGrad(Y), vTurbulenceGrad(Z), "
           << "vDirection(X), vDirection(Y), vDirection(Z), "
           << "Magnitude, "
           << "vTurbPQR(P), vTurbPQR(Q), vTurbPQR(R), " << endl;
    }
    if (from == 2) {
      cout << vTurbulenceNED << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
    }
  }
  if (debug_lvl & 64) {
    if (from == 0) { // Constructor
      cout << IdSrc << endl;
      cout << IdHdr << endl;
    }
  }
}



} // namespace JSBSim