ramanterr(386).cpp

处理激光雷达数据

/*********************************************************/
/* PROG for calculating air density and Temperature (k)  */
/*  with N2 Raman Lidar                                  */
/*  Hn and Hc are not same                               */
/*  Many groups of raw MCS data are ok (0005101.386)     */
/*                                                       */
/*  Ver.3.0,  Iternational calculating for density       */
/*                                                       */
/* Edited on 2000-7-30                                   */
/*  gz[k]=9.8*6371.0/(6371.0+high[k])*6371.0/(6371.0+high[k]) */
/* Corrected on 2001-01-08                               */
/* to calculate statistic errors                         */
/* On 2001-07-28                                         */
/*********************************************************/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<process.h>
#define HH 10         /* smooth signal, higher and lower height-areas */
#define Deltazh 0.9   /* h> HH km */
#define Deltazl 0.45   /* h<hh km */

#define HH2 10         /* smooth result  */
#define HH3 15         /* smooth result  */
#define Deltazh3 1.0   /* h> HH3 km */
#define Deltazh2 0.60   /* h> HH2 km */
#define Deltazl2 0.3   /* h<HH2 km */

#define  MCS  -4
#define  B1 850   /* start channelNo for calcuating background */
#define  B2 990   /* end channelNo for calcuating background */

/************************************************/
/* Main Prog for  air density and temperature   */
/***********************************************/
void main()
{
  void readmcs(char *, float *, float *,int *,float *);
  FILE *fp1,*fp2,*fp3,*fsig;
  FILE *faer,*fo3,*ftau,*fact;
  int i,j,k,nm,nn,gatechan,times;
  float layer,mypb,allpb;  // layer= km
  int   newnum,pointc,pointnc,pointtc;
  float tmph[1000],high[1000], hp[1000],p[1000],p2[1000]; // para. for raw data
  float mh[610],mdensity[610],mt[610],mp[610]; /* para. for air model data */
  float ah[610],aerex[610],r[610];  /* para. for aerosol (L300 data) */
   float am[600],ap[600];
  float taumol,tauaer,taum[600],taua[600],tauo[600],tauall[600];
  float x[600],density[600],snr[600],nerr[600],column[600],tp[600],terr[600];
                                  // para. for calculating air density
  float density2[600],tp2[600];
  float snrc,w[40],ww,allpp,deltaz;    // para. for smoothing
  float Hnc,Nc, Htc,Ntc,Tc,PZ2c,tmp,aven,flag1,flag2;
  char name[18],st;
  float gz[600], avegz, wv, TIMES,denlidar, den0;
  float gefh[600],gefact[600];
  layer=0; newnum=0; mypb=0; allpb=0; gatechan=0; tmp=0;
  denlidar=0; den0=0;
  wv=1.2;        /* aerosol wavelength exponent between 0.5- 2.0 */
  TIMES=2;       /* iteration order: 1, 2, 5  */

  Hnc=3.0;     /* Normalized height for air density  */
  flag1=1;

  Htc=5;   /* Reference height for temperature, where SNR>=5  */
  flag2=1;

  for(i=0;i<1000;i++)
     { tmph[i]=0; high[i]=0; p[i]=0; hp[i]=0;  snr[i]=0; }
  for(i=0;i<610;i++)
    { mh[i]=0; mt[i]=0; mp[i]=0; mdensity[i]=0;  }
  for(i=0;i<600;i++)
    { x[i]=0;  snr[i]=0; density[i]=0;  column[i]=0;
      tp[i]=0; nerr[i]=0; terr[i]=0;
      taum[i]=0;  taua[i]=0;  tauo[i]=0;  tauall[i]=0;
     am[i]=0;   ap[i]=0;  density2[i]=0; tp2[i]=0;
      gz[i]=0.0; gefh[i]=0;gefact[i]=0;
    }

 if( (fp1=fopen("atmod355.dat","r"))==NULL ) /* atmos.model with N(z),T(z) */
    { printf("\n Cannot open mol355.dat \n");
      exit(1);
    }
 if( (fp2=fopen("tresult.dat","w"))==NULL)    /* Temperature result file */
    { printf("\n Cannot open lidart.dat\n");
       exit(1);
    }
 if( (fp3=fopen("lidarn.dat","w"))==NULL)   /* Air density result file */
    { printf("\n Cannot open lidarn.dat\n");
       exit(1);
    }
 if( (fsig=fopen("sig.dat","w"))==NULL) /* Signal result file */
    { printf("\n cannot open sig.dat\n");
      exit(1);
    }
 if( (faer=fopen("aer355_98.dat","r"))==NULL) /* Aerosol data file */
    { printf("\n Cannot open aer.dat \n");
      exit(1);
    }

 if( (ftau=fopen("ftau.dat","w"))==NULL) /* Output single-trip trans. */
    { printf("\n Cannot open ftau.dat \n");
      exit(1);
    }
  if((fact=fopen("fact.dat","r"))==NULL)
  { printf("\n Cannot open fact.dat \n");
     exit(1);
  }
for(; ;)
   {
    printf("\n (1) Input high-altitude filename: ");
    scanf("%18s",name);
    fflush(stdin);

    readmcs(name,p,&layer,&newnum,&mypb); /* P[i] already subtract Pb */
    allpb=allpb+mypb;
    for(i=0;i<newnum;i++)
      { tmph[i]=(i+1)*layer;    /* layer unit: km */
        hp[i]=p[i]+hp[i]; /* to intergrate all of groups */
       // printf("%d  %6.2f %6.0f\n", i, high[i], hp[i]);
      }
    printf("\n (2) Enter another high datafile:( Y/N ?) ");
    fflush(stdin);
    st=(char)getchar();
    if(st=='n'||st=='N') break;
   }   /* end reading high-altitude file */

  /* to skip gate or overlapping height */
  k=0;
  printf("(3) Input start channelNo(gate): ");
  scanf("%d", &gatechan);
  fflush(stdin);

 /* to skip gatechannel height or overlapping height */
 for(i=0; i<newnum;i++)
    { if(i>=gatechan)
         { if(hp[i]<=0 )
             { printf(" *** high=%6.2f  P=%6.2f\n",tmph[i],hp[i]);
               break;
             }
           p2[k]=hp[i];
           high[k]=tmph[i];
           snr[k]=hp[i]/sqrt(double(hp[i]+allpb));
      if( fabs(high[k]-4)<=0.015) printf("H=%6.2f  SNR=%6.2f\n",high[k],snr[k]);
      if( fabs(high[k]-6)<=0.015) printf("H=%6.2f  SNR=%6.2f\n",high[k],snr[k]);
      if( fabs(high[k]-8)<=0.015) printf("H=%6.2f  SNR=%6.2f\n",high[k],snr[k]);
      k++;
        }
    }
 newnum=k;   /* new channel number */

 /* to output signal with subtracting background noise */
 fprintf(fsig,"H(km)    Pz-Pb     P(z)*z*z \n");
 for(i=0; i<newnum; i++)
   fprintf(fsig,"%6.2f  %6.2f   %6.3e\n", high[i],p2[i],p2[i]*high[i]*high[i]);
 fclose(fsig);
 if(newnum>=600) newnum=600;

//* mean-running smoothing  raw signal */
 k=0;
 for(i=0;i<newnum;i++)
   {
     if(high[i]<=HH ) nn=Deltazl/layer/2; /* lower areas */
     else nn=Deltazh/layer/2;            /* higher areas */
     if(i>=nn && i<newnum-nn)
       { allpp=0;
	      for(j=-nn;j<=nn;j++)
     	      allpp=allpp+p2[i+j]/(2*nn+1);
         p[k]=allpp;
         high[k]=high[i];
         snr[k]=snr[i]*sqrt(2*nn+1);
         k++;
       }
	 else
	 {
		 p[k]=p2[i];
		 k++;
	 }
	 }
newnum=k;
for(i=0;i<600;i++)
{
	fscanf(fact,"%f %f\n",&gefh[i],&gefact[i]);
}
 for(i=0;i<newnum;i++)
   { for(j=0;j<600;j++)
       { if(high[i]>=gefh[j]&& high[i]<gefh[j+1])
           {
	     gefact[i]=gefact[j]+( gefact[j+1]-gefact[j])/(gefh[j+1]-gefh[j])*(high[i]-gefh[j]);
          
            break;
           }
        }
 }
for(i=0;i<newnum;i++)
{
  if(high[i]<4.0)
 p[i]=p[i]/gefact[i];
}


for(i=0;i<newnum;i++)
  {
    if( fabs(high[i]-Hnc)<=0.015) /* to find more exact calibrate height */
        { pointnc=i;
          Hnc=high[i];
        }
     if( fabs(high[i]-Htc)<=0.015) /* to find more exact calibrate height */
       { pointtc=i;
         Htc=high[i];
         break;
       }
   }
printf(" NewCali.:Hnc=%6.2f  Htc=%6.2f  %d  %d %d  HNC=%6.3f HTC=%6.3f\n",
         Hnc,Htc,newnum,pointnc,pointtc,high[pointnc],high[pointtc]);

 /***********************************************************/
 /* to read air model or NMC temperarure and peressure data */
 /***********************************************************/
 /* to read model N2 and T(K) data, ME90TN2.dat */
 for(i=70;i>=0;i--)
  {
   fscanf(fp1,"%f %f %f %f \n",&mh[i],&mt[i],&mp[i],&mdensity[i]);
    density[i]=7.247249e18*mp[i]/mt[i]; 
   //  printf("mh=%6.3f  %6.3e\n",mh[i],mdensity[i]);
  }     /* mh[i]/km,  mdensity[i]/molcule/cm3  */
 fclose(fp1);

 /*To find Reference temp(K) at Htc(km) */
for(j=0;j<70;j++)        /* Note:  difference between Halo and NMC/Sage data */
  {
   if(Htc>=mh[j] && Htc<mh[j+1])    /* for HALOE or CIRA86 model data */
	   { Tc=mt[j]+(mt[j+1]-mt[j])/(mh[j+1]-mh[j])*(Htc-mh[j]);
        printf("Htc=%6.3f   Tc=%6.4e\n",Htc,Tc);
        break;
     }
  }

/* interpolate mol. density data at Lidar height (km) */
 for(i=0;i<newnum;i++)
   { for(j=0;j<70;j++)
       { if(high[i]>=mh[j]&& high[i]<mh[j+1])
           {am[i]=density[j]+( density[j+1]-density[j])/(mh[j+1]-mh[j])*(high[i]-mh[j]);
           // am[i]=am[i]*(2.776701e-26)*(1.0e+5);  /* molex./km */
            break;
           }
        }
    // printf("Mol H=%6.3f  am=%6.4e \n",high[i],am[i]);
    }
 printf (" OK 1  (Model density N(z), T(K) and Mol.ex,  %d\n",newnum);

/***************************************************/
/* to get molecule,aerosol    */
/***************************************************/
/* Aerosol extinction */
/* to read aerosol data: aer355.dat 98 year */
 for(i=0;i<610;i++)
   fscanf(faer,"%f %f %f\n",&ah[i],&r[i],&aerex[i]);
       /* ah[i]/km,  aerex[i]/ km-1  */
 fclose(faer);
 /* interpolate aerosol data at Lidar height (km) */
 for(i=0;i<newnum;i++)
   { for(j=0;j<600;j++)
       { if(high[i]>=ah[j]&& high[i]<ah[j+1])
           {ap[i]=aerex[j]+( aerex[j+1]-aerex[j])/(ah[j+1]-ah[j])*(high[i]-ah[j]);
            //ap[i]=0.0;  /* ignore aerosol influence */
            break;
           }
        }
    // printf("Aer H=%6.3f %6.4e \n",high[i],am[i]);
   }
printf (" OK 2 (aerosol.ext) %d\n",newnum);

/*******************************************/
/* to calculte  single-trip transmission   */
/******************************************/
times=0;
for(; ;)
{
 /*To find Reference density  at Hnc(km) */
 for(j=0;j<newnum;j++)       /* Note:  difference between Halo and NMC/Sage data */
  {
   if(Hnc>=high[j]&&Hnc<high[j+1])
     { Nc=am[j]+(am[j+1]-am[j])/(high[j+1]-high[j])*(Hnc-high[j]);
       printf("Hnc=%6.3f   Nc=%6.4e\n",Hnc,Nc); /* Nc: /cm3 */
       break;
     }
  }

den0=0;
for(j=0;j<newnum;j++)  /* calcuate mol.ex */
 { den0=den0+am[j];    /* column density */
   am[j]=am[j]*(2.776701e-26)*(1.0e+5);  /* molex./km */
  }
printf("Times: %d  den0=%6.6e\n",times+1,den0);

printf(" Printnc=%d  Printtc=%d  Hnc=%6.3f  h=%6.3f  layer=%6.3f \n",
           pointnc,pointtc,Hnc,high[pointnc],layer);
fprintf(ftau,"H(km)   am     ap       taum     taua       tauall\n");
for(i=0;i<newnum;i++)
  { taumol=0; tauaer=0; 
    if(high[i]<Hnc)
      { for(j=i;j<=pointnc;j++)
          { taumol=taumol+am[j]*(high[j+1]-high[j]);
            tauaer=tauaer+ap[j]*(high[j+1]-high[j]);
          }
        taum[i]=(taumol-(am[i]+am[pointnc])/2*layer)*(1+(float)pow(355.0/386,4.0));
        taua[i]=(tauaer-(ap[i]+ap[pointnc])/2*layer)*(1+(float)pow(355.0/386,wv));
     
        taum[i]=exp(-taum[i]);  /* molecule trans. <=1.0 */
        taua[i]=exp(-taua[i]);  /* aerosol trans. <=1.0 */
      }
    else if( fabs(high[i]-Hnc)<=0.015)
           { taum[i]=1.0, taua[i]=1.0, tauo[i]=1.0; }
    else
       {for(j=pointnc;j<=i;j++)
          { taumol=taumol+am[j]*(high[j+1]-high[j]);
            tauaer=tauaer+ap[j]*(high[j+1]-high[j]);
         
          }
        taum[i]=(taumol-(am[i]+am[pointnc])/2*layer)*(1+(float)pow(355.0/386,4.0));
        taua[i]=(tauaer-(ap[i]+ap[pointnc])/2*layer)*(1+(float)pow(355.0/386,1.2));