ra532_607b.cpp

处理激光雷达数据

/***************************************************************************/
/* PROG for calculating aerosol scattering ratio                           */
/*  1. N2 Raman and Mie signals                                            */
/*  2. Input *.532 and *.607 MCS data                                      */
/*                                                                         */
/* Edited on 2004-9-10                                                    */
/* No Channel geometrical factor correction                                                                        */
/***************************************************************************/

#include "readmcs_1.c"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <process.h>
#include <string.h>



void main()
{
  FILE *fpout,*fmol,*faer, *fr;
   int i,j,k;

  char fname[20];

 float ratio1[1000];

 for(i=0;i<1000;i++)
  {
    ratio1[i]=0;
   }


/********************************************************************/
/*  Read MCS file  "*.532","*.607"
/*  Read channel geometrical correction MCS file, 532channel and
/*       607-channel.  calculate factor ,
/*   And make correction
/********************************************************************/

 float p532[1000],p607[1000], temp[1000];
  char str1[20],str2[20],str3[20],filename[20];

 for(i=0;i<1000;i++)
  {
   p532[i]=0; p607[i]=0; temp[i]=0;
   }

 printf("Input date,such as: 040315\n");
 scanf("%s",str1);

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

 printf("\n\nInput MCS file (.532) number,such as: 2 or 3 or 4....\n");
 scanf("%d",&k);


 for(i=1;i<k+1;i++)
  {
     strcpy(filename,"");;
     itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"00");
     strcat(filename,str2);
     strcat(filename,".532");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"%s\n",filename);
     for(j=0;j<1000;j++)
      p532[j]=p532[j]+temp[j];
   }

  /* for(j=0;j<1000;j++)
     printf("%d  %.3f\n",j,p532[j]);   */

  printf("\n\n Input MCS file (n?.607) number,such as: 2 or 3 or 4....\n");
  scanf("%d",&k);


  for(i=1;i<k+1;i++)
  {
     strcpy(filename,"");;
     itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"00");
     strcat(filename,str2);
     strcat(filename,".607");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"\n%s",filename);
     for(j=0;j<1000;j++)
        p607[j]=p607[j]+temp[j];
   }

 /*for(j=0;j<1000;j++)
   fprintf(fpout,"\n%d  %.3f",j,p607[j]);         */


 fprintf(fpout,"\n\n");
 for(j=0;j<1000;j++)
   fprintf(fpout,"%d  %.4e    %.4e\n",j,p532[j],p607[j]);
fclose(fpout);

/*********************************************************
**
**   subtract background and data smoothing
**
**********************************************************/
    int start_bg, end_bg, gate, nn, newnum;
  float ave532=0, ave607=0,  spa_res, pp1, pp2, layer;
  float p1[1000],p2[1000],high[1000];



  printf("\n Input bakground range from channel ? to channel ?, such as: 500  700\n");
  scanf("%d%d",&start_bg,&end_bg);
  fflush(stdin);

  for(i=start_bg;i<end_bg;i++)
     {
      ave532=ave532+p532[i]/(end_bg-start_bg);
      ave607=ave607+p607[i]/(end_bg-start_bg);
      }

  k=0;
  printf("\n Input start channel No(gate): ?  and spatial resolution(m),such as:150 m \n");
  scanf("%d%f", &gate,&spa_res);      /* gatechannel start channel */
  layer=spa_res/1000.;
  fflush(stdin);
  for(i=0;i<1000;i++)
     {
       if(i>=gate)
          {  //if(p532[i]<=0 || p607[i]<=0) break;
            high[k]=i*spa_res/1000.;
            p532[k]=p532[i]-ave532;
            p607[k]=p607[i]-ave607;
 //           fprintf(fsig,"%6.3f    %6.3f    %6.3f\n",high[k],p532[k],p607[k]);
            k++;
         }
     }
  newnum=k;

/* to smooth lidar signal */
 for(i=0;i<newnum;i++)
    {
      p1[i]=p532[i]; p2[i]=p607[i];
      if(high[i]<2.5 ) nn=0.3/layer/2; /* lower areas */
      else nn=0.9/layer/2.0;            /* higher areas */

      if(i>=nn && i<(newnum-nn))
        { pp1=0;  pp2=0;
	       for(j=-nn;j<=nn;j++)
    	      { pp1=pp1+p532[i+j]/(2*nn+1);
              pp2=pp2+p607[i+j]/(2*nn+1);
             }
         p1[i]=pp1;  p2[i]=pp2;
       }
     if( p1[i]<=0 || p2[i]<0) break;
    }
 newnum=i;

 for(i=0;i<newnum;i++)
   {
    ratio1[i]=p1[i]/p2[i];
    p532[i]=p1[i];
    p607[i]=p2[i];
    }



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

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

 if( (fmol=fopen("atmod532r.dat","r"))==NULL ) // atmos.model with N(z),T(z)
    { printf("\n Cannot open atmod532r.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[1000];

 for(i=0;i<newnum;i++)
   { for(j=0;j<nm-1;j++)
       { 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;
     }
  newnum=i;


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

   float taum[1000],scale[1000], scale0=0, H1, H2, tmp;
   for(i=0;i<1000;i++)
   {
	   taum[i]=0;scale[i]=0;
   }

  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*(high[i]-high[i-1]);
    tmp=taum[i]*(1-pow(532.0/607,4.0));
    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]);
  }

/* to find out scale-height */
float Hc, p532c, p607c;
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];
           p532c=p532[i];    p607c=p607[i];
          }
       }
    if(high[i]>H2) break;
  }
printf("Hc=%6.2f  pointc=%d  p532c=%6.2f  p607c=%6.2f  scale0=%6.2f\n",
          Hc, pointc, p532c, p607c, scale0);
printf("SNR532(Hc)=%6.2f   SNR607(Hc)=%6.2f\n",sqrt(p532c),sqrt(p607c));

/* 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];
  for(i=0;i<500;i++)
  { ah[i]=0;
    aerex[i]=0;
	ap[i]=0;}
   /* 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<100;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;
           }
        }
   }

/*******************************************/
/* to calculte  scattering ratio R(z)     */
/******************************************/
  float r[1000],molt[1000],aert[1000];
  r[pointc]=1.0; molt[pointc]=1.0;  aert[pointc]=1.0;
  float tmp1=0,  tmp2=0, tmp3=0;
  float deltaz, aveam, aveap=0, newap;
  float Cbh, Cth, v, s1, Rc, pHmax;
  pHmax=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(532.0/607.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=50; v=1.0; }

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

      r[j]=ratio1[j]*p607c/p532c*Rc* molt[j]* aert[j];
      newap=(r[j]-1)*am[j]*3/8/3.14*s1;
      //printf("newap=%6.4e\n",newap);

      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=50; v=1.0; }
           tmp3=tmp3+deltaz*aveap* (1-(float)pow(532.0/607.0,v));
           aert[j]=(float)exp(-tmp3);
           r[j]=ratio1[j]*p607c/p532c*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(532.0/607.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=50; v=1.0; }
     tmp2=tmp2-deltaz*aveap*(1-(float)pow(532.0/607.0,v));
     aert[j]=(float)exp(-tmp2);

     r[j]=ratio1[j]*p607c/p532c*Rc*molt[j]* aert[j];
     newap=(r[j]-1)*am[j]*3/8/3.14*s1;

     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=50; v=1.0; }
           tmp3=tmp3-deltaz*aveap* (1-(float)pow(532.0/607.0,v));
           aert[j]=(float)exp(-tmp3);

           r[j]=ratio1[j]*p607c/p532c*Rc* molt[j]* aert[j];
           newap=(r[j]-1)*am[j]*3/8/3.14*s1;
           tmp2=tmp3;
           k++;
         }
   }

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

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

 fprintf(fr,"H/km   R(532)    Bp/km\n");
 for(i=0; i<newnum; i++)
   { tmp1=r[i];
     if(high[i]>=1.5)
      { if(high[i]>=1.5 && high[i]<=3.0) nn=0.30/layer/2;
        else if(high[i]>3.0 && high[i]<=6.0) nn=0.6/layer/2;
        else nn=0.6/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;
    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);

 }