386_355a.cpp

处理激光雷达数据

/***************************************************************************/
/* PROG for calculating aerosol scattering ratio                           */
/*  1. N2 Raman and Mie signals                                            */
/*  2. Input *.355 and *.386 MCS data                                      */
/*                                                                         */
/* Edited on 2004-03-29                                                    */
/*   by yanshunsheng  2007.9.24                                                                      */
/***************************************************************************/

#include "READMCS_1.C"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <process.h>
#include <string.h>



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

  char fname[20];

 float ratio1[1000],ratio2[1000];

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


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

 float p355[1000],p386[1000],pnb[1000],pnc[1000], temp[1000];
  char str1[20],str2[20],str3[20],filename[20];

 for(i=0;i<1000;i++)
  {
   p355[i]=0; p386[i]=0;pnb[i]=1;pnc[i]=1; temp[i]=0;
 // printf("%f %f %f %f %f\n", p355[i],p386[i],pnb[i],pnc[i], temp[i]);
 }

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

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

 printf("\n\nInput MCS file (.355) 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,".355");
     printf("%s   ",filename);

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

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

  printf("\n\n Input MCS file (n?.386) 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,"n00");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

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

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

 printf("\n\n Channel geometrical factor correction?\n");
 printf("     Yes, choise  '1' \n");
 printf("      No, choise  '0' \n");
 scanf("%d",&p);

 if(p==1)
  {
   for(i=0;i<1000;i++)
    {
      pnb[i]=0;
      pnc[i]=0;
     }
   printf("\n Input channel geometric MCS file (c?.386) number\n");
   printf("       such as: 1 or 2 or 3... \n");
   scanf("%d",&k);

    fprintf(fpout,"\n");
   for(i=1;i<k+1;i++)
    {
     strcpy(filename,"");;
     itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"wc200");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

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

   printf("\n Input channel geometric MCS file (b?.386) number\n");
   printf("       such as: 1 or 2 or 3... \n");
   scanf("%d",&k);

   fprintf(fpout,"\n");
   for(i=1;i<k+1;i++)
    {
     strcpy(filename,"");;
      itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"wb200");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"\n%s",filename);
     for(j=0;j<1000;j++)
        pnb[j]=pnb[j]+temp[j];
     }
 // for(j=0;j<1000;j++)
    //  fprintf(fpout,"%d  %.3f\n",j,pnb[j]);    
  }

 fprintf(fpout,"\n\n");
 for(j=0;j<1000;j++)
   fprintf(fpout,"%f %.4e    %.4e    %.4e     %.4e\n",j*0.03,p355[j],p386[j],pnc[j],pnb[j]);
fclose(fpout);
//printf("%d  %.4e    %.4e    %.4e     %.4e\n",j*0.03,p355[j],p386[j],pnc[j],pnb[j]);


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



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

  for(i=start_bg;i<end_bg;i++)
     {
      ave355=ave355+p355[i]/(end_bg-start_bg);
      ave386=ave386+p386[i]/(end_bg-start_bg);
      avenc=avenc+pnc[i]/(end_bg-start_bg);
      avenb=avenb+pnb[i]/(end_bg-start_bg);
      }

  k=0;
  printf("\n Input start channel No(gate): ?  and spatial resolution(m),such as:30m \n");
  scanf("%d%f", &gate,&spa_res);      /* gatechannel start channel */
  layer=spa_res/1000.;// layer=0.03 km
  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.;
            p355[k]=p355[i]-ave355;
            p386[k]=p386[i]-ave386;
            pnc[k]=pnc[i]-avenc;
            pnb[k]=pnb[i]-avenb;
            p355[i]=p355[k];//
            p386[i]=p386[k];
            pnc[i]=pnc[k];
            pnb[i]=pnb[k];
            high[i]=high[k];
			//printf("%6.3f    %6.3f    %6.3f\n",high[i],p355[i],p386[i]);//两杠间2007.3.20改过
		   k++;
         }
     }
  newnum=k;

/* to smooth lidar signal */
 for(i=0;i<newnum;i++)
    {
      p1[i]=p355[i]; p2[i]=p386[i];//2007.3.20日p355[i]改为 p355[k]
      if(high[i]<2.5 ) nn=0.3/layer/2; /* lower areas */
      else nn=0.6/layer/2.0;            /* higher areas */

      if(i>=nn && i<(newnum-nn))
        { pp1=0;  pp2=0;
	       for(j=-nn;j<=nn;j++)
    	      { pp1=pp1+p355[i+j]/(2*nn+1);// 300 m or 600 m smooth
              pp2=pp2+p386[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++)
   {
    p355[i]=p1[i];
    p386[i]=p2[i];
	 
	 ratio1[i]=p1[i]/p2[i];
    

    printf("%f \n",ratio1[i]);
 }
   
   
/* to smooth channel sigal */

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

 for(i=0;i<newnum;i++)
    {
      p1[i]=pnc[i]; p2[i]=pnb[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+pnc[i+j]/(2*nn+1);
            pp2=pp2+pnb[i+j]/(2*nn+1);
            }
        p1[i]=pp1;  p2[i]=pp2;
      }
     if( p1[i]<=0 || p2[i]<0) 
		 break;
   }

 for(i=0;i<newnum;i++)
   if(high[i]<10)
     ratio2[i]=p1[i]/p2[i];

/***********************************************************************/

 float ave_cal=0, cal_high;

 printf("Input calibration range, such as: 2 or 3 or 4\n");
 scanf("%f",&cal_high);

 j=0;
 for(i=0;i<newnum;i++)
    if(high[i] > cal_high && high[i] < cal_high+1)
     {
       ave_cal=ave_cal+ratio2[i];
       j++  ;
      }
 ave_cal=ave_cal/j;

 strcpy(str3,"");
 strcat(str3,str1);
 strcat(str3,"fact_aer.dat");

 if( (fact=fopen(str3,"w"))==NULL)
    {
     printf("\n cannot open %s\n",str3);
     exit(1);
    }
 for(i=0;i<newnum;i++)
   if(high[i]<10)
     {
      ratio2[i]=ratio2[i]/ave_cal;
      fprintf(fact,"%.3f   %.3f\n", high[i],ratio2[i]);
      }