co2386_371a.cpp

处理激光雷达数据

/***************************************************************************/
/* PROG for calculating co2  mixing ratio                           */
/*  1. N2 and co2 Raman  signals                                         */
/*  2. Input *.371 and *.386 MCS data                                      */
/*   using 98 extinction  mode                                                                      */
/* Edited on 2004-03-31                                                    */
/*                                                                         */
/***************************************************************************/

#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,*fraw,*fs,*fw;
   int i,j,k,p;

//  char fname[20];

 float ratio1[1000],ratio2[1000];

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


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

 float p371[1000],p386[1000],pnb[1000],pna[1000], temp[1000];
  char str1[20],str2[20],str3[20],str4[20],str5[20],str6[20],filename[20];

 for(i=0;i<1000;i++)
  {
   p371[i]=0; p386[i]=0;pnb[i]=1;pna[i]=1; temp[i]=0; 
   }

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

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

 printf("\n\nInput MCS file (w?.407) group 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,"co200");
     strcat(filename,str2);
     strcat(filename,".371");
     printf("%s   ",filename);

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

  /* for(j=0;j<1000;j++)
     printf("%d  %.3f\n",j,p407[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,"l00");
     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,p386[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;
      pna[i]=0;
     }
   printf("\n Input channel geometric MCS file (co2b11?.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,"co2b1100");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

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

   printf("\n Input channel geometric MCS file (co2c12?.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,"co2c1200");
     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,"%d  %.4e    %.4e    %.4e     %.4e\n",j,p371[j],p386[j],pna[j],pnb[j]);
fclose(fpout);

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

 strcpy(str5,"");
 strcat(str5,str1);
 strcat(str5,"平滑原始信号out.dat");

 if( (fs=fopen(str5,"w"))==NULL)
    {
     printf("\n cannot open %s\n",str5);
     exit(1);
    }
fprintf(fs,"高度    二氧化碳拉曼光子数   氮气拉曼光子数   标定二通道       三通道\n  ");
  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++)
     {
      ave371=ave371+p371[i]/(end_bg-start_bg);
      ave386=ave386+p386[i]/(end_bg-start_bg);
      avena=avena+pna[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:30 \n");
  scanf("%d%f", &gate,&spa_res);      /* gate channel start channel */
  layer=spa_res/1000.;
  fflush(stdin);
  for(i=0;i<1000;i++)
     {
       if(i>=gate)
          {  //if(p407[i]<=0 || p386[i]<=0) break;
            high[k]=i*spa_res/1000.;
            p371[k]=p371[i]-ave371;
            p386[k]=p386[i]-ave386;
            pna[k]=pna[i]-avena;
            pnb[k]=pnb[i]-avenb;
          //fprintf(fs,"%6.3f    %6.3f    %6.3f     %f     %f\n",high[k],p371[k],p386[k],pna[k],pnb[k]);
            k++;
         }
     }
  newnum=k;

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

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

 strcpy(str6,"");
 strcpy(str6,str1);
 strcat(str6,"未被修正outco2.dat");
 if( (fw=fopen(str6,"w"))==NULL)
    {
     printf("\n cannot open %s\n",str3);
     exit(1);
    }

 for(i=0;i<newnum;i++)
   { 
	 if(high[i]<10)
    ratio1[i]=p1[i]/p2[i];
    p371[i]=p1[i];
    p386[i]=p2[i];
   fprintf(fw,"%.5f   %.5f\n", high[i],ratio1[i]);

    } 

fclose(fw);
/* to smooth channel sigal */

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

 for(i=0;i<newnum;i++)
    {
      p1[i]=pna[i]; p2[i]=pnb[i];
      if(high[i]<2.5 ) nn=0.15/layer/2; /* lower areas */
      else nn=0.3/layer/2.0;            /* higher areas */

      if(i>=nn && i<(newnum-nn))
        { pp1=0;  pp2=0;
	       for(j=-nn;j<=nn;j++)
    	      { pp1=pp1+pna[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)
     {
  fprintf(fs,"%6.3f    %6.3f    %6.3f     %f     %f\n",high[i],p371[i],p386[i],pna[i],pnb[i]);
      }
  fclose(fs);
 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;
  printf("%f\n",ave_cal);

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

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


 for(i=0;i<newnum;i++)
 {
	 if(high[i]<3)
 ratio1[i]=ratio1[i]/ratio2[i];
 }
 strcpy(str4,"");
 strcpy(str4,str1);
 strcat(str4,"co2raw_result.dat");
 if( (fraw=fopen(str4,"w"))==NULL)
    {
     printf("\n cannot open %s\n",str3);
     exit(1);
    }
    
 for(i=0;i<newnum;i++)
   if(high[i]<10)
     {
      
      fprintf(fraw,"%.4f   %.4f\n", high[i],ratio1[i]);
      }
  fclose(fraw);

//*************************************************************************
//
//   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[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]);
          
           
        }
      if (high[i]>=molHmax) break;
     }
  newnum=i;

 /* read 1998 mean aerosol data     */

                                                           // 首次猜测气溶胶文件名
  if( (faer=fopen("aer355_98.dat","r"))==NULL) /* Aerosol data file */
    { printf("\n Cannot open file aer355_98.dat \n");
      exit(1);                                                                   // !!!!
    }


  float ah[1000],aerex[1000],ap[1000];
   /* Aerosol extinction */

  i=0;
  while( (fscanf(faer,"%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<330;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]);
           
      }

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

   float taum[1000],scale[1000], tmp1, tmp2, taua[1000];

  for(i=0;i<1000;i++)
  {
	  taum[i]=0;
    scale[i]=0;
     taua[i]=0;
  }
  for(i=1;i<newnum;i++)
  {
    taum[i]=taum[i-1]+(am[i]+am[i-1])/2.0*(high[i]-high[i-1]);
    tmp1=taum[i]*(1-pow(354.7/386,4.0))-taum[i]*(1-pow(354.7/371,4.0));

    taua[i]=taua[i-1]+(ap[i]+ap[i-1])/2.0*(high[i]-high[i-1]);
    tmp2=taua[i]*(1-354.7/386)-taua[i]*(1-354.7/371);

    scale[i]=ratio1[i]*exp(-tmp1-tmp2);
  }

/*******************************************/
/* to calculte co2 vapor mixing ratio   */
/*******************************************/
  float cons;

  printf(" \n If calibrated by radiosonde, choise '1'.\n");
  printf(" else  with a constant, choise '2'.\n");
  scanf("%d",&k);

  if(k==1)
     printf("ok\n");
    else
     {
      printf("Input calibration constant:?\n");
      scanf("%f",&cons);
      }

  for(i=1;i<newnum;i++)
     {
      scale[i]=scale[i]/cons;
     }

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

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

 fprintf(fr,"Height/km   co2_vapor_mixing_ratio/g/kg\n");
 for(i=1; i<newnum; i++)
   { tmp1=scale[i];
     if(high[i]>=1.5)
      { if(high[i]>=1.5 && high[i]<=3.0) nn=0.150/layer/2;
        else if(high[i]>3.0 && high[i]<=6.0) nn=0.3/layer/2;
        else nn=0.3/layer/2;
        if(i>=nn&& i<newnum-nn)
      	 { tmp1=0;
	         for(j=-nn;j<=nn;j++)
      	      tmp1=tmp1+scale[i+j]/(2*nn+1);
           }
      }
    scale[i]=tmp1;

    if(high[i]<=13)
      fprintf(fr,"%4.3f   %6.5f\n",high[i],scale[i]);
  }

 fclose(fr);
 printf("\n ! co2 vapor mixing ratio: %s\n", filename);


 }