?? 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);
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -