ramanterr(386).cpp

处理激光雷达数据

    

        taum[i]=exp(taum[i]);
        taua[i]=exp(taua[i]);
      
       }
   tauall[i]= taum[i]*taua[i];
   fprintf(ftau,"%6.3f  %6.4e  %6.4e  %6.4e   %6.4e   %6.4e  \n",
             high[i],am[i],ap[i],taum[i],taua[i],tauall[i]);
  }

flag1=1;
/* 4. To  compute air density...... */
for(i=0;i<newnum;i++)
   { x[i]=p[i]*high[i]*high[i];       /* to caculate range-corrected signals */
     if(flag1==1)
       if( fabs(high[i]-Hnc)<=0.015)  /* normalizing height for air densty */
           { PZ2c=x[i];
             flag1=-1;
           }
    }
printf(" Hnc=%6.3f  Nc=%6.4e   PZ2C=%6.4e \n",Hnc,Nc,PZ2c);

fprintf(fp3,"H(km)     N(/cm3)       G(z)\n");
for(i=0;i<newnum;i++)
   {
    density[i]=x[i]/PZ2c*Nc*tauall[i];     /* molecule density: unit:  /cm3 */
    gz[i]=9.8*6371.0/(6371.0+high[i])*6371.0/(6371.0+high[i]);

    nerr[i]=1.0/snr[i];
    printf("%6.3f  %6.6e\n",high[i],nerr[i]);
   }
/* to smooth air density, due to lower SNR */ 
for(i=0; i<newnum;i++)
  {
   if(high[i]<=HH2 ) nn=Deltazl2/layer/2; /* lower areas */
   else nn=Deltazh2/layer/2;            /* higher areas */

   density2[i]=density[i];
   if(i>=nn && i<(newnum-nn))
      { tmp=0;
        for(j=-nn;j<=nn;j++)
   	     tmp=tmp+density[i+j]/(2*nn+1);
        density2[i]=density[i];
	    }
    fprintf(fp3,"%3.3f  %6.5e  %6.3f\n",high[i],density2[i],gz[i]);
  }

 denlidar=0;
 for(i=0;i<newnum;i++)
    denlidar=denlidar+density2[i]; /* unit: /cm3,Lidar density column */

 //printf("Times: %d  denlidar=%6.6e  Deln=%6.6f\n",times+1,denlidar,(denlidar-den0)/den0);
 times=times+1;
 //printf("\n Iteration times: %d\n",times);
 if( times>=TIMES) break;
 else
   {
    for(i=0;i<newnum;i++)
      am[i]=density2[i];  /* new iteration  */
   }
 }  /* end iteration density */
 fclose(ftau);
 fclose(fp3);

 flag2=1;
 for(i=0;i<newnum;i++)
   { density2[i]=density2[i]*(4.65040718667e-20);  /* unit: Kg/m3, Shibata method */
     if(flag2==1)  /* to get air density for temperature */
       if( fabs(high[i]-Htc)<=0.015)
         { Ntc=density2[i];
           pointc=i;
           flag2=-1;
         }
   }

 /* 5. to  compute air temperature  with Shibata method  */
 for( i=0;i<newnum;i++)
   { tmp=0;
     if(high[i]<Htc)
      { for(j=i;j<pointc;j++)
         { deltaz=(high[j+1]-high[j]);    /* unit:Kilometer */
           aven=(density2[j]+density2[j+1])/2.0;
           avegz=(gz[j]+gz[j+1])/2.0;

           tmp=deltaz*aven*avegz;// positive value when z < zc */
           column[i]=column[i]+tmp;
         }
       }
     else if(high[i]==Htc) column[i]=0;
     else
      { for(j=pointc;j<i;j++)
         { deltaz=(high[j+1]-high[j]);    /* unit:Kilometer */
           aven=(density2[j]+density2[j+1])/2.0;
           avegz=(gz[j]+gz[j+1])/2.0;

           tmp=-deltaz*aven*avegz;   /* negative value when Z> Zc*/
           column[i]=column[i]+tmp;
          }
      }
    tp[i]=(Tc*Ntc+28.966/8.31*column[i])/density2[i]; /* calculating temperature */

    terr[i]=(tp[i]- 28.966/8.31*gz[i]/2)* nerr[i];   /* unit: K */
   }

 k=0;
 for(i=0; i<newnum ;i++)
   {
    if(high[i]>=0.6)    // to choose 7<High<25 km 'result
       { high[k]=high[i];       snr[k]=snr[i];
         density2[k]=density2[i]; nerr[k]=nerr[i];
         tp[k]=tp[i];            terr[k]=terr[i];
         k++;
        }
    }
 newnum=k;

 fprintf(fp2,"H(km)   T(K)   terr\n");
 for(i=0; i<newnum;i++)
  { if(high[i]<=HH2 ) nn=Deltazl2/layer/2.0; /* lower areas */
    else if( high[i]>=HH2 && high[i]<=HH3) nn=Deltazh2/layer/2.0; /* Mid areas */
    else  nn=Deltazh3/layer/2.0;  /* highest areas */

    tp2[i]=tp[i];
    if(i>=nn && i<newnum-nn)
       { tmp=0;
         for(j=-nn;j<=nn;j++)
   	     tmp=tmp+tp[i+j]/(2*nn+1);
         tp2[i]=tmp;
	    }
    if(high[i]<25.0)
     fprintf(fp2,"%3.2f   %6.3f    %6.2f\n",high[i],tp2[i],terr[i]);
   }
 fclose(fp2);
 printf("\n ** Output signal file: sig.dat.\n");
 printf("\n ! Temperature result file: lidart.dat.\n");
 printf("\n ! Air density result file: lidarn.dat.\n");

 }  /* to end all the program */


/* function readmcs() to read the raw data */
/*****************************************************************/
/*        Sample program compatible with Microsoft and Borland C */
/*	   to read header and channel data from a MCS data file 	     */
/* to output format: channel No.- Photon counter                 */
/*****************************************************************/
void readmcs(char *pname, float *pp,float *player,int *point,float *pbb)
 {
  char 	  reserved[100];  /* buffer to skip reserved bytes */
 short int		  f_type;	  /* MCS file type,2 bytes,Windows 98/Pen3 computer*/
   char 	  trigger;	  /* Trigger Flag */
   char		  dwell_flag;	  /* External Dwell Flag */
   char		  dwell_units;	  /* 0=us, 1=ms, 2=sec, 3=ns */
   char		  acq_mode;	  /* Acquisition mode flag 0=replace, 1=sum */
   unsigned long  dwell_913;	  /* Dwell time in old 913 format */
 unsigned short int   pass_length;    /* pass length in channels */
   unsigned long  pass_count;
   unsigned long  pass_count_preset;
   char 	  acq_time[9];    /* buffer for acquisition time */
   char 	  acq_date[9];    /* buffer for acquisition date */
 unsigned short int   mark_chan;	  /* first marker channel */
   char		  mcs_num;        /* 1-8 are valid  */
   char		  cal_flag;	  /* 0=no calibration */
   char		  cal_units[4];	  /* calibration units */
   float	  cal_zero;	  /* calibration zero intercept */
   float	  cal_slope;	  /* calibration slope */
   char		  id_byte;	  /* always 0xaa */
   char		  detector_len;	  /* Detector description length */
   char		  detector[64];	  /* detector description */
   char		  sample_len;	  /* Sample description length */
   char		  sample[64];	  /* Sample description */
   char		  disc_sel;	  /* 0=SCA otherwise Disc */
   char		  disc_edge;	  /* 0=falling else rising */
   float	  disc;		  /* disc setting in volts */
   float	  sca_uld;	  /* sca upper-level in volts */
   float	  sca_lld;	  /* sca lower-level in volts */
   float	  dwell;	  /* dwell time in dwell_units */
   char		  consistant;	  /* settings consistant flag */
   char		  mcs_id[9];	  /* MCS ID string "0914-001" */
   unsigned long  chan_data;      /* stores channel data */

unsigned short int   count;          /* general purpose loop counter  */
   FILE  *f_pointer;

   int k;

   if( (f_pointer = fopen(pname,"rb"))==NULL)
      { printf(" Cannot open file %s", pname);
        exit(1);
      }

/***************************************************/
/*           Header Data                           */
/*		    Read header info from MCS file			   */
/***************************************************/
/* Read filetype -4 (MCS) */

   fread(&f_type,sizeof(short int),1,f_pointer);
   if (f_type != MCS) {
      printf("Not a valid file\n");
      exit(1);
   }
   fread(&trigger,sizeof(char),1,f_pointer);	/* Read Trigger Flag */
   fread(&dwell_flag,sizeof(char),1,f_pointer);	/* Read dwell flag */
   fread(&dwell_units,sizeof(char),1,f_pointer);/* Read dwell units */
   fread(&acq_mode,sizeof(char),1,f_pointer);
   fread(&dwell_913,sizeof(long),1,f_pointer);

   fread(&pass_length,sizeof(unsigned short int),1,f_pointer);
   fread(&pass_count,sizeof(long),1,f_pointer);
   fread(&pass_count_preset,sizeof(long),1,f_pointer);
   fread(acq_time,sizeof(char),8,f_pointer);
   fread(acq_date,sizeof(char),8,f_pointer);
   fread(&mark_chan,sizeof(unsigned short int),1,f_pointer);
   fread(&mcs_num,sizeof(char),1,f_pointer);
   fread(&cal_flag,sizeof(char),1,f_pointer);
   fread(cal_units,sizeof(char),4,f_pointer);
   fread(&cal_zero,sizeof(float),1,f_pointer);
   fread(&cal_slope,sizeof(float),1,f_pointer);
   fread(reserved,sizeof(char),10,f_pointer);
   fread(&id_byte,sizeof(char),1,f_pointer);
   fread(reserved,sizeof(char),1,f_pointer);
   fread(&detector_len,sizeof(char),1,f_pointer);
   fread(detector,sizeof(char),63,f_pointer);
   fread(&sample_len,sizeof(char),1,f_pointer);
   fread(sample,sizeof(char),63,f_pointer);
   fread(reserved,sizeof(char),16,f_pointer); /* skip view info & reserved */
   fread(&disc_sel,sizeof(char),1,f_pointer);
   fread(&disc_edge,sizeof(char),1,f_pointer);
   fread(&disc,sizeof(float),1,f_pointer);
   fread(&sca_uld,sizeof(float),1,f_pointer);
   fread(&sca_lld,sizeof(float),1,f_pointer);
   fread(&dwell,sizeof(float),1,f_pointer);
   fread(&consistant,sizeof(char),1,f_pointer);
   fread(reserved,sizeof(char),21,f_pointer);
   fread(mcs_id,sizeof(char),8,f_pointer);

   mcs_id[8]=0;
   printf("MCS%2i,  MCS_ID = %8s\n",mcs_num+1,mcs_id);
   printf("Pass Length= %5i, Pass Count = %4i, Pass Count Preset = %4i\n",
       pass_length,pass_count,pass_count_preset);

   if (dwell_flag == 0)
    {
      printf("Dwell Time = %2.2g",dwell);
      if (dwell_units==0)
        { printf("us\n");  *player=dwell*300/2/1000;} /* unit:km */
      else if (dwell_units==1)
        { printf("ms\n"); *player=dwell*300000/2/1000; }
      else if (dwell_units==2)
        { printf("s\n");  *player=dwell*300000000/2/1000; }
      else
        { printf("ns\n"); *player=dwell*0.3/2/1000; }
    } else printf("Dwell Time = External\n");

   if (trigger==0) printf("Trigger = Internal, ");
   else printf("Trigger = External, ");

   if (acq_mode==0) printf(" Acquisition Mode = Replace\n");
   else printf(" Acquisition Mode = Sum\n");

   acq_time[8]=0;
   acq_date[8]=0;
   printf("Data collected at %8s on %8s\n",acq_time,acq_date);

 /* if (cal_flag==0) printf("Spectrum is not Calibrated\n");
   else printf("Calibration Slope = %13.7g,  Offset = %13.7g,  units=%.4s\n",
               	cal_slope,cal_zero,cal_units);

   if (detector_len==0) printf("No hardware description\n");
   else {
      printf("Hardware Description:\n");
      for (count=0; count < (unsigned)detector_len; count++)
        putchar(detector[count]);
      printf("\n");
   }

   if (sample_len==0) printf("No sample description\n");
   else {
      printf("Sample Description:\n");
      for (count=0; count < (unsigned)sample_len; count++)
         putchar(sample[count]);
      printf("\n");
     }

   if (disc_sel == 0) {
      printf("SCA Selected:  LLD=%13.7gV,  ULD=%13.7gV\n",sca_lld,sca_uld);
      }
   else {
      printf("Disc Selected: Voltage=%13.7gV,  ",disc);
      if (disc_edge==0) printf("Edge=Falling\n");
      else printf("Edge=Rising\n");
     }
   if (consistant==0)
      printf("WARNING:  Settings are not consistant with data\n"); /*

/************************************************/
/*             Channel Data                     */
/*		    Output channel data from MCS file     */
/************************************************/
  *pbb=0; k=0;
  for (count = 0; count < pass_length; count++)
    {
      fread(&chan_data,sizeof(long),1,f_pointer);
      *(pp+count)=chan_data;
      if(count>B1 && count<B2) /* to calculte background noise */
        { *pbb= (*pbb)+chan_data;
          k++;
         }
    }
  fclose(f_pointer);

  *pbb=(*pbb)/k;  /* averaged value as background noise */
  printf("Pbb=%6.4f \n",*pbb);
  for (count = 0; count < pass_length; count++)
     *(pp+count)= (*(pp+count))- (*pbb);
  *point=pass_length;
   }  /* to end readmcs() */