aerallerr(低中高信号).cpp

处理激光雷达数据

   }
 ww=x2[k]/bm[k];
 kk=0;
 for(i=k; high[i]<h2-layer; i++)
   { if(x2[i]/bm[i]<=ww)
   	{ ww=x2[i]/bm[i];
	     kk=i;
	     h1=high[i];
      }
    }

 if(choice==1) r[kk]=1.01;   /* at 532 nm */
 if(choice==2) r[kk]=1.002;  /* at 355 nm */

 b[kk]=(r[kk]-1.0)*bm[kk];  // aer'backscatter coffecient.
 printf("\n * Hc=%4.2fkm  SNR=%4.1f  S1=%3.1f Rc=%4.2f*\n",h1,snr[kk],s1,r[kk]) ;

 /* 12. Calculate aerosol scattering ratio from order to order */
 for( i=kk;i>=1;i--)  // forward to interger
    { aa=(s1-s2)*(bm[i]+bm[i-1])*layer;
      mm=x[i-1]*exp(aa);
      b[i-1]=mm/(x[i]/(b[i]+bm[i])+s1*(x[i]+x[i-1]*exp(aa) )*layer)-bm[i-1];
      r[i-1]=1.0 +b[i-1]/bm[i-1];
     //printf("-- %3.1f  %6.2e  %6.2e  %6.2e\n",high[i],r[i],b[i],b[i]*s1);
    }
 for( i=kk;i<newnum-1;i++)   // backward to interger
    { aa=(s1-s2)*(bm[i]+bm[i+1])*layer;
      mm=x[i+1]*exp(-aa);
      b[i+1]=mm/(x[i]/(b[i]+bm[i])-s1*(x[i]+x[i+1]*exp(-aa) )*layer)-bm[i+1];
      r[i+1]=1.0+b[i+1]/bm[i+1];
    //  printf("++ %3.1f  %6.2f  %6.2e  %6.2e\n",high[i],r[i],b[i],b[i]*s1);
    }

 /* 13. smooth and output results (0.45-1.0 km)  */
 fprintf(fout,"H(km)   SR    Aerex(/km)    rerr%    bperr%\n");
 for(i=0; i<newnum-1; i++)
   { pp1=r[i];
     pp2=b[i];
     if(high[i]>=5)
      { if(high[i]>=5 && high[i]<=20) nn=0.6/layer/2;   // 450m for smoothing
        else nn=1.05/layer/2;
        if(i>=nn&& i<newnum-nn)
      	 { pp1=0; pp2=0.0;
	         for(j=-nn;j<=nn;j++)
      	     { pp1=pp1+r[i+j]/(2*nn+1);
   	           pp2=pp2+b[i+j]/(2*nn+1);
	            }
  	        }
      }
    r[i]=pp1;
    b[i]=pp2;

    bperr[i]=rerr[i]*r[i]/(r[i]-1);

   if( high[i]<=35.0)
     fprintf(fout,"%4.2f  %6.3f  %6.3e  %6.2f  %6.2f\n",high[i],r[i],b[i]*s1,
                    rerr[i]*100,bperr[i]*100);

  }
 fclose(fp2);    fclose(fout);

 printf("\n Signal file: L625p.dat.\n");
 if(choice==1)
    printf("\n\n ! Result file: err532.dat.\n");
 else if(choice==2)
   printf("\n\n ! Result file: a355.dat.\n");
 else
    printf("\n\n! Error result file.\n");

  }  /* to main prog */



/*****************************************************************/
/*        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)
 {
  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;

   float pb;
   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     */
/************************************************/
  pb=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 */
        { pb=pb+chan_data;
          k++;
         }
    }
  fclose(f_pointer);

  pb=pb/k;  /* average value for background noise */
  printf("Pb=%6.5f \n",pb);
  for (count = 0; count < pass_length; count++)
     *(pp+count)= (*(pp+count))-pb;

  *point=pass_length;

 }  /* to end readmcs() */