aerallerr(低中高信号).cpp

处理激光雷达数据

/***************************************************************/
/* To calculate aerosol'scattering ratio and extinction coffi. */
/*  with Lidar by Fernald methold in Appl.Opt 23,1984,P652     */
/*  1. At 532 nm and 355 nm   (S1=50 @532 nm, S1=30@ 355nm     */
/*  2. Using low and high-altitude OR only high-data           */
/*       (to junct a whole signals)                            */
/*  3. Many groups of data are ok                              */
/*  4. To calculate statistic error SQRT(p(z))                 */

/*  Edited on 2000-3-22, tested by WU, Status: ok              */
/*  Corrected on 2001-6-24                                     */
/***************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<process.h>

#define  MCS  -4
#define  B1 700  /* start channelNo for calcuating background */
#define  B2 950  /* end channelNo for calcuating background */

void main()
 {
   void readmcs(char *, float *, float *,int *);
   FILE *fout, *fp2, *fsig;
   char name[16],st;
   float tmph[1000],high[1000],p[1000],hp[1000],lp[1000],mp[1000],p2[1000],snr[1000];
   float layer,hpp,lpp,mpp,coeff,coeff1 ;
   int i,k,newnum, lowchan, midchan1,midchan2,highchan, gatechan;

   int j,choice;
   int   nn,w[20];  /* para for smoothing signals */
   float ww,allpp; /* para for smoothing signals */
   int nm,kk;         /* para for reading molecule data */
   float hm[100],am[100],tm[100],tmol[100],taum, s1,s2;  /* molecule data */
   float r[270],b[270],x[270],bm[270],x2[270], rerr[270],bperr[270];
   float  h1,h2, aa, mm,pp1,pp2, molhmax; /* scale height range */

   for(i=0;i<100;i++)
    { hm[i]=0; am[i]=0; tm[i]=0; tmol[i]=0; }
   for(i=0;i<270;i++)
     { bm[i]=0;  x[i]=0.0;   r[i]=0;  b[i]=0; rerr[i]=0; bperr[i]=0; }

   layer=0; newnum=0;
   for(i=0;i<1000;i++)
     { tmph[i]=0; high[i]=0; p[i]=0; hp[i]=0; lp[i]=0; mp[i]=0;p2[i]=0; snr[i]=0; }
   if( (fsig=fopen("L625p.dat","w"))==NULL)
	   { printf("\n cannot open L625P.dat\n");
	     exit(1); /*  mod-atmosph,mol-scatter&extinction */
	   }
 for(; ;)
   {
    printf("\n (1) Input high-altitude filename: ");
    scanf("%14s",name);
    fflush(stdin);
    readmcs(name,p,&layer,&newnum);

    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 */

  printf("\n (3) You use mid-altitude datafile:( Y/N ?) ");
  fflush(stdin);
  st=(char) getchar();
  if(st=='Y'||st=='y')
   { for(; ;)
     {
       printf("\n (3) Input mid-altitude filename: ");
       scanf("%14s",name);
       fflush(stdin);
       readmcs(name,p,&layer,&newnum);

       for(i=0;i<newnum;i++) /* to intergrate all of groups data */
         { tmph[i]=(i+1)*layer;
           mp[i]=p[i]+mp[i]; /* to intergrate all of groups */
         }
       printf("\n (4) Enter another mid datafile:( Y/N ?) ");
       fflush(stdin);
       st=(char) getchar();
       if(st=='n'||st=='N') break;
       } /* end reading mid-altitude file */
  printf("\n (3) You use low-altitude datafile:( Y/N ?) ");
  fflush(stdin);
  st=(char) getchar();
  if(st=='Y'||st=='y')
   { for(; ;)
     {
       printf("\n (3) Input low-altitude filename: ");
       scanf("%14s",name);
       fflush(stdin);
       readmcs(name,p,&layer,&newnum);

       for(i=0;i<newnum;i++) /* to intergrate all of groups data */
         { tmph[i]=(i+1)*layer;
           lp[i]=p[i]+lp[i]; /* to intergrate all of groups */
         }
       printf("\n (4) Enter another low datafile:( Y/N ?) ");
       fflush(stdin);
       st=(char) getchar();
       if(st=='n'||st=='N') break;
       }
      mpp=0; lpp=0;
      printf("(5.1) Input overlap channelNo(low  mid1):  ");
      scanf("%d %d", &lowchan, &midchan1);
      fflush(stdin);
      for(i=0; i<newnum; i++) /* to calculate  */
       { if(i>=lowchan && i<=midchan1)
         { mpp=mpp+mp[i];  lpp=lpp+lp[i]; }
        }
      coeff=mpp/lpp;
      printf(" coeff=%6.4f\n", coeff);
   } 
  else { lowchan=-1; highchan=-1; }
     mpp=0; hpp=0;
     printf("(5.2) Input overlap channelNo(mid2 high): ");
	 scanf("%d  %d", &midchan2, &highchan);
      fflush(stdin);
      for(i=0; i<newnum; i++) /* to calculate  */
       { if(i>=midchan2 && i<=highchan)
         { mpp=mpp+mp[i];  hpp=hpp+hp[i]; }
        }
      coeff1=mpp/hpp;
      printf(" coeff1=%6.4f\n", coeff1);
   } /* end reading low-altitude datafile */
  else { lowchan=-1; highchan=-1; }
    
  k=0;
  printf("(6) Input start channelNo(gate): ");
  scanf("%d", &gatechan);
  fflush(stdin);
 /* to get a whole signal cureve */
 for(i=0; i<newnum;i++)
    { if(i>=gatechan)
       { if(i<lowchan)
         {p2[k]=lp[i]*coeff;  high[k]=tmph[i];  }
         else if (i>=lowchan && i<=midchan1)
		 { p2[k]=(lp[i]*coeff+mp[i])/2;  high[k]=tmph[i];}
		 else if (i>midchan1 && i<=midchan2)
		 {  p2[k]=mp[i]; high[k]=tmph[i];}
		  else if (i>midchan2 && i<=highchan)
		  { p2[k]=(mp[i]+hp[i]*coeff1)/2;  high[k]=tmph[i];
		  }
         else
            { if(hp[i]<=0 ) break;
              p2[k]=hp[i]*coeff1; high[k]=tmph[i];
            }
         snr[k]=(float)sqrt(double(p2[k]));

        if( fabs(high[k]-25)<0.1) printf("H=%6.2f  SNR=%6.1f\n",high[k],snr[k]);
        if( fabs(high[k]-30)<0.1) printf("H=%6.2f  SNR=%6.1f\n",high[k],snr[k]);
        k++;
        }
    }
 newnum=k;   /* new channel number */
 if(newnum>=270) newnum=270;

 /* output signal */
 for(i=0; i<newnum; i++)
     fprintf(fsig,"%6.2f  %6.2f\n", high[i],p2[i]);
 fclose(fsig);

 printf("\n  *************************************\n");
 printf("  * (7) Wavelength choice:            *\n");
 printf("  *                        1-532nm    *\n");
 printf("  *                        2-355nm    *\n");
 printf("  *************************************\n");
 scanf("%d",&choice);   fflush(stdin);
  if(choice==1)
     {
		if( (fp2=fopen("atmod532.dat","r"))==NULL)
		   { printf("\n cannot open atmod532.dat\n");
		     exit(1); /*  mod-atmosph,mol-scatter&extinction */
		   }
		if( (fout=fopen("err532.dat","w"))==NULL)
		    { printf("\n cannot open err532.dat\n");
		      exit(1);  /* result file for r[] ,aerbs[],molbs[] */
		    }
      s1=50;  /* ext to backscatter ratio of are. ?40 */
		}
  else if(choice==2)
     {
	   if( (fp2=fopen("atmod355.dat","r"))==NULL)
		   { printf("\n cannot open atmod355.dat\n");
		      exit(1); /* mod-atmosph,mol-scatter&extinction */
		   }
		if( (fout=fopen("a355.dat","w"))==NULL)
		    { printf("\n cannot open a355.dat\n");
		       exit(1);  /* result file for r[] ,aerbs[],molbs[] */
		     }
      s1=30;       /* ext to backscatter ratio of are. ?40 */
    }
 else { printf("Error choice.\n"); exit(1);  }

 /* 8. to smooth raw photon-counters */
 nn=0;
 printf("\n  (9) Enter the number of smoothing point: ");
 scanf("%d",&nn);  /* the number of data points, odd: 1 3 5 3 1 */
 fflush(stdin);
 for( i=0;i<nn/2; i++)
   { w[i]=2*i+1;
     w[nn-i-1]=2*i+1;  /* triangle-weight coffe. */
     p[i]=p2[i];  /* without smoothing */
   }
 for(i=nn/2;i<newnum-nn/2;i++)
   { ww=0;  allpp=0;
     for(j=0;j<nn;j++)
       { allpp=allpp+p2[i-nn/2+j]*w[j];
         ww=ww+w[j];
       }
     p[i]=allpp/ww;
     rerr[i]=1.0/snr[i]/sqrt(nn);
   }

 /* 9. read data of "atmod532.dat",hm[i],am[] */
 s2=8*3.1415926/3;        /* ext to backscatter ratio of molecule */
 fscanf(fp2,"%d",&nm);
 for(i=0;i<nm;i++)
   fscanf(fp2,"%f %f %f %f %f\n",&hm[i],&am[i],&taum,&taum,&taum);
      /* am[i] mol'ext. at 532nm */
 molhmax=hm[i-1];
 tm[0]=0;

 for(i=1;i<nm;i++)
   tm[i]=tm[i-1]+(am[i]+am[i-1])/2.0*(hm[i]-hm[i-1]);

 /* printf("\n To begin to join mol and aer'ext and scat........\n"); */
 for(i=0;i<newnum; i++)
  { for(j=0;j<nm;j++) /* to interpolate mol'scatter & ext coeffi. */
      if(high[i]>=hm[j] && high[i]<hm[j+1])
   	 {
	     bm[i]=am[j]+( am[j+1]-am[j])/(hm[j+1]-hm[j]) *(high[i]-hm[j]);
        bm[i]=bm[i]/s2;           /* 532 nm */
        tmol[i]=tm[j]+( tm[j+1]-tm[j])/(hm[j+1]-hm[j])*(high[i]-hm[j]);
       }
    if(high[i]>=molhmax) break;
  }
newnum=i;

 /* 10. to begin to compute aer'r[],ab[],mb[]......... */
 for(i=0;i<newnum;i++)
    {
     x[i]=p[i]*high[i]*high[i];
     x2[i]=x[i]*(float)exp(2*tmol[i]);
     // fprintf(fout,"%4.2f %6.4f\n",high[i],1.0/snr[i]);
    }

 printf("\n  (11) Enter scale-height range (h1 h2): ");
 scanf("%f %f",&h1,&h2);  fflush(stdin);
 for(i=0;i<newnum;i++)
   { if(high[i]>=h1 && high[i]<h2)
       { k=i; break; }