386_355a.cpp

处理激光雷达数据

  fclose(fact);

 for(i=0;i<newnum;i++)
 { if(high[i]<cal_high)
	 ratio1[i]=ratio1[i]/ratio2[i];
  printf("%f \n",  ratio1[i]);//标定高度下修正ratio1
 }
	



//*************************************************************************
//
//   read atmospheric model  "atmod355.dat".
//
//************************************************************************

   int nm;
 float mh[500], molex[500], mdensity, mp, mt, molHmax;

 if( (fmol=fopen("atmod355.dat","r"))==NULL ) // atmos.model with N(z),T(z)
    { printf("\n Cannot open atmod355.dat \n");
      exit(1);
    }

 fscanf(fmol,"%d\n", &nm); /* to read model N2 and T(K) data, ME90TN2.dat */
 for(i=0;i<nm;i++)
  {                                                              //读取氮气消光,密度
   fscanf(fmol,"%f %e %f %f %f\n",&mh[i],&molex[i],&mdensity,&mp,&mt);
  }     /* mh[i]/km,  molex[i]/ km-1  */
 molHmax=mh[i-1];
 fclose(fmol);

/* interpolate molecular data at Lidar height (km) */
 float am[1500];

 for(i=0;i<newnum;i++)
   { for(j=0;j<nm;j++)//原为j<nm,会导致错误
       { if(high[i]>=mh[j]&& high[i]<mh[j+1])
           {am[i]=molex[j]+( molex[j+1]-molex[j])/(mh[j+1]-mh[j])*(high[i]-mh[j]);
           
        break;
           }
        }
      if (high[i]>=molHmax) break;
     printf("high=%6.2f     am= %6.3f  \n",   high[i], am[i]);      //error
     }
  newnum=i;

 
/******************************************************************************/
/* to calculate molcular optical depth (0-z) */

  float taum[1000]={0} ,scale[1000], scale0, H1, H2, tmp;

  printf("\nPlease enter Cali.height(H1, H2 km):  ");
  scanf("%f %f",&H1,&H2);
  printf("H1=%6.1f    H2=%6.1f\n",H1,H2);                    //标高H1, H2
 //H1=6.0;   /* Down-boundary height km */
 // H2=12.0;  /* Up-boundary height km */

 
  for(i=1;i<newnum;i++)
  {
    taum[i]=taum[i-1]+(am[i]+am[i-1])/2.0*layer;//layer=0.03,积分
    tmp=taum[i]*(1-pow(355.0/386.0,4));
    scale[i]=ratio1[i]*exp(-tmp);
    if(fabs(high[i]-H1)<=layer) scale0=scale[i];
    printf("high=%6.2f  taum=%6.4f  scale=%6.3f\n",high[i],taum[i],scale[i]);
    printf("ratio=%6.2f  \n",ratio1[i]);
  }

/* to find out scale-height */
float Hc, p355c, p386c;
int pointc;
printf("scale0=%6.4f\n\n",scale0);
for(i=0;i<newnum;i++)
  { if(high[i]>=H1 && high[i]<=H2)
      { if(scale[i]<=scale0)
          {scale0=scale[i];  pointc=i;  Hc=high[i];
           p355c=p355[i];    p386c=p386[i];
          }
       }
    if(high[i]>H2) break;
  }
printf("Hc=%6.2f  pointc=%d  p355c=%6.2f  p386c=%6.2f  scale0=%6.2f\n",
          Hc, pointc, p355c, p386c, scale0);
printf("SNR355(Hc)=%6.2f   SNR386(Hc)=%6.2f\n",sqrt(p355c),sqrt(p386c));

/* read aerosol data     */

 // printf("Please enter first-guess aerosol file:  ");
 // scanf("%s",fname);                                          // 首次猜测气溶胶文件名
//  if( (faer=fopen(fname,"r"))==NULL) /* Aerosol data file */
  //  { printf("\n Cannot open file %s \n",fname);
  //    exit(1);                                                // !!!!
  //  }

 // char ttmp[20];
 // float ah[500],aerex[500],ap[500];
   /* Aerosol extinction */
 // fscanf(faer,"%s %s %s\n",ttmp,ttmp,ttmp);
//  i=0;
//  while( (fscanf(faer,"%f %f %f\n",&ah[i],&tmp,&aerex[i]))!=EOF)
 //     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<nm-1;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]);
             
		//	break;
         //  }
      //  }
 // printf("high=%6.2f     ap= %6.3f  \n",   high[i], ap[i]);  //error?
 //}

 printf("Please enter first-guess aerosol file:  ");
  scanf("%s",fname);                                          // 首次猜测气溶胶文件名
  if( (faer=fopen(fname,"r"))==NULL) /* Aerosol data file */
    { printf("\n Cannot open file %s \n",fname);
      exit(1);                                                // !!!!
    }

  char ttmp[20];
  float ah[500],aerex[500],ap[500];
   /* Aerosol extinction */
  fscanf(faer,"%s %s \n",ttmp,ttmp);
  i=0;
  while( (fscanf(faer,"%f %f %f\n",&ah[i],&aerex[i]))!=EOF)
      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<1000;j++)//原为j<nm, aer.dat只能计算1.95km 以后的消光系数,可能导致ratio1在1.95前不正常,平滑更导致全体错误
       { 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]);
            break;
           }
        
  }
   printf("high=%6.2f     ap= %6.3f  \n",   high[i], ap[i]);
 }
/*******************************************/
/* to calculte  scattering ratio R(z)     */
/******************************************/
 float r[1000],molt[1000],aert[1000];
 float hh[1000],rr[1000],apr[1000];
  r[pointc]=1.0; molt[pointc]=1.0;  aert[pointc]=1.0;
  float tmp1=0,  tmp2=0, tmp3=0;
  float deltaz, aveam, aveap, newap;
  float Cbh, Cth, v, s1, Rc, pHmax;
  int a;
  pHmax=18;//18
  Cbh=0.0;    /* Cloud base height */
  Cth=0.0;  /* Cloud top height */
  Rc=1.002;//可改

 for(j=(pointc-1);j>=0;j--)     /* High< Hc */
    { deltaz=(high[j+1]-high[j]);      /* unit:Kilometer */
      aveam=(am[j+1]+ am[j])/2.0;
      tmp1=tmp1+deltaz*aveam* (1-(float)pow(355.0/386.0,4.0));
      molt[j]=(float)exp(-tmp1);

      aveap=(ap[j+1]+ ap[j])/2.0;

      if(high[j]>=Cbh && high[j]<=Cth) { s1=10; v=0.0; }
      else { s1=30; v=1.0; }//原为30

      tmp2=tmp2+deltaz*aveap* (1-(float)pow(355.0/386.0,v));
      aert[j]=(float)exp(-tmp2);

      r[j]=ratio1[j]*p386c/p355c*Rc* molt[j]* aert[j];
      newap=(r[j]-1)*am[j]*3/8/3.14*s1;
    //printf("newap=%6.4e\n",newap);
 printf("j=%d %6.3f  %6.4f %6.4f %6e %6.3f\n",j,high[j],r[j],ratio1[j],aert[j],molt[j]);//aert[j]:error
      k=0;
      while( fabs((newap-ap[j])/newap)>0.01)
         { ap[j]=newap;
           aveap=(ap[j+1]+ ap[j])/2.0;

           if(high[j]>=Cbh && high[j]<=Cth) { s1=10; v=0.0; }
           else { s1=30; v=1.0; }
           tmp3=tmp3+deltaz*aveap* (1-(float)pow(355.0/386.0,v));
           aert[j]=(float)exp(-tmp3);
           r[j]=ratio1[j]*p386c/p355c*Rc* molt[j]* aert[j];

           newap=(r[j]-1)*am[j]*3/8/3.14*s1;
          // printf("**** newap=%6.4e\n",newap);
           tmp2=tmp3;
           k++;
         }
   //printf("j=%d %6.3f  %6.4e  %6.4e\n",j,high[j],tmp1,molt[j] );
     }
 tmp1=0;   tmp2=0;  tmp3=0;
 for(j=(pointc+1);j<newnum;j++)   // H>=HC
   { deltaz=(high[j]-high[j-1]);      /* unit:Kilometer */
     aveam=(am[j]+ am[j-1])/2.0;
     tmp1= tmp1 -deltaz*aveam*(1-(float)pow(355.0/386.0,4.0));
     molt[j]=(float)exp(-tmp1);

     aveap=(ap[j]+ ap[j-1])/2.0;

     if(high[j]>=Cbh && high[j]<=Cth) { s1=10; v=0.0; }
     else { s1=30; v=1.0; }
     tmp2=tmp2-deltaz*aveap*(1-(float)pow(355.0/386.0,v));//tmp2-deltaz中间为什么变为减号
     aert[j]=(float)exp(-tmp2);

     r[j]=ratio1[j]*p386c/p355c*Rc*molt[j]* aert[j];
     

     newap=(r[j]-1)*am[j]*3/8/3.14*s1;
	printf("J=%d %6.3f  %6.4f %6.4f %6e %6.3f\n",j,high[j],r[j],ratio1[j],aert[j],molt[j]);
    
	 k=0;
     while( fabs((newap-ap[j])/newap)>0.01)
        {  ap[j]=newap;
           aveap=(ap[j]+ ap[j-1])/2.0;
 
           if(high[j]>=Cbh && high[j]<=Cth) { s1=10; v=0.0; }
            else { s1=30; v=1.0; }
            tmp3=tmp3-deltaz*aveap* (1-(float)pow(355.0/386.0,v));
            aert[j]=(float)exp(-tmp3);

            r[j]=ratio1[j]*p386c/p355c*Rc* molt[j]* aert[j];
            newap=(r[j]-1)*am[j]*3/8/3.14*s1;//newap为消光系数
            tmp2=tmp3;
            k++;
          }
    }

 /*  smooth and output results (0.45-1.0 km)  */

 strcpy(filename,"");
 strcat(filename,str1);
 strcat(filename,"r355.dat");
 if( (fr=fopen(filename,"w"))==NULL)
    {
     printf("\n cannot open %s\n", filename);
     exit(1);
    }

 fprintf(fr,"H/km   R(355)    Bp/km\n");
 //for(i=0; i<newnum; i++)
//if(high[i]<=pHmax&& high[i]>=1.95)
     // fprintf(fr,"%4.2f   %6.4f     %6.4e\n",high[i],r[i],ap[i]);
//for(i=newnum; i>=0; i--)

//a=0;
//while( (fscanf(fr,"%f %f %f\n",&hh[a],&rr[a]))!=EOF)
 //     a++;   

//{
//	for( a=0;a<newnum; a++)

//	tmp1=rr[a];
     
 //      { if(hh[a]>=1.95 && hh[a]<=3.0) nn=0.30/layer/2;// 300 m or 600 m smooth
 //      else if(hh[a]>3.0 && hh[a]<=6.0) nn=0.6/layer/2;
 //       else nn=0.6/layer/2;
 //      if(a>=nn&& a<newnum-nn)
 //    	 { tmp1=0;
//	        for(j=-nn;j<=nn;j++)
  //    	      tmp1=tmp1+rr[a+j]/(2*nn+1);
 //         }
 //     }
 //   rr[a]=tmp1;
   //apr[s]=(rr[s]-1)*am[i]*3/8/3.14;
 //   if(hh[a]<=pHmax)
 //     fprintf(fr,"%4.2f   %6.4f   \n",hh[a],rr[a]);
 // }
for(i=0; i<newnum; i++)
{    tmp1=r[i];
    if(high[i]>=1.98)	
	{  // if (high[i]=1.95)  nn=0;
		 if(high[i]>=1.98&& high[i]<2.07)  nn=1;
		 if(high[i]>=2.10 && high[i]<=3.0) nn=0.30/layer/2;//nn=5
	     
       else if(high[i]>3.0 && high[i]<=18) nn=0.6/layer/2;
        else nn=1;//.5/layer/2;
        if(i>=nn&& i<=newnum-nn)
      	 { tmp1=0;
	         for(j=-nn;j<=nn;j++)
     	      tmp1=tmp1+r[i+j]/(2*nn+1);
          }
     }
  r[i]=tmp1;
   ap[i]=(r[i]-1)*am[i]*3/8/3.14;//此中ap为气溶胶后向散射系数
    if(high[i]<=pHmax)
      fprintf(fr,"%4.2f   %6.4f     %6.4e\n",high[i],r[i],ap[i]);
 }


 fclose(fr);
 printf("\n ! aerosol scattering ratio: %s\n", filename);
 }