利用CE318太阳光度计数据处理反演AOD和水汽含量wv

#include <stdio.h> #include <math.h> //#define M 179 #define TMLOGPATH "E:\\DTF5\\CE318data_software\\x081111.txt" #define V0RSDPATH "E:\\DTF5\\CE318data_software\\V0_R_SD.txt" #define DEPTHPATH "E:\\DTF5\\CE318data_software\\AODx081111.txt" #define ALLPATH "E:\\DTF5\\CE318data_software\\allx081111.txt" #define ROW 17 void main() { int i,j,M=0; FILE *p; char ch,str[10]="x081111"; double timal[ROW],m[ROW],V[ROW][5],mavg,a[5],b[5],lat=31.9037,lon=117.1584;//大营盘lat=40.255,lon=115.8;//合肥lat=31.9037,lon=117.1584;风电场lat=41.017,lon=114.758//经纬度注意是否要修改 double Vavg[5],V0[5],total=0.0,sumV[5],dx=0.0,dy[5],dxy[5],R[5],SD[5]; double bm[ROW],bmavg,Vm[ROW],Vmavg; double tao[ROW][5],taor[5],taoa8[ROW][5],taoa940[ROW],tao940[ROW],taoH2O[ROW],dxy940=0.0,ab=0.606,ba=0.580;//,Wb[M],W[M] double taoa5[ROW][4],lnnamta5[4],namta5avg,sumtaoa5[ROW],taoa5avg[ROW],dxym[ROW],slope1[ROW],arfa[ROW],beta[ROW]; double namta8[5]={1.020,0.936,0.870,0.670,0.440}; for(i=0;i<5;i++) { dxy[i]=0.0; dy[i]=0.0; sumV[i]=0.0; SD[i]=0.0; } for(i=0;i<ROW;i++) { dxym[i]=0.0; sumtaoa5[i]=0.0; } if((p=fopen(TMLOGPATH,"r"))==NULL) { printf("can't open the file!\n"); return; } ch=fgetc(p); while(ch!=EOF) { if(ch=='\n') M++; ch=fgetc(p); } fclose(p); printf("%d\n",M); if((p=fopen(TMLOGPATH,"r"))==NULL) { printf("can't open the file!\n"); return; } for(i=0;i<M;i++) { fscanf(p,"%lf %lf",&timal[i],&m[i]); for(j=0;j<5;j++) fscanf(p,"%lf",&V[i][j]); } fclose(p); //计算大气层外仪器测量值V0及定标因子C for(i=0;i<M;i++) total=total+m[i]; mavg=total/M; //printf("%lf\n",mavg); for(i=0;i<M;i++) for(j=0;j<5;j++) sumV[j]=sumV[j]+V[i][j]; for(j=0;j<5;j++) { Vavg[j]=sumV[j]/M; //printf("%lf\t",exp(Vavg[i])); } for(i=0;i<M;i++) dx=dx+pow((m[i]-mavg),2); for(i=0;i<5;i++) for(j=0;j<M;j++) { dy[i]=dy[i]+pow((V[j][i]-Vavg[i]),2); dxy[i]=dxy[i]+(m[j]-mavg)*(V[j][i]-Vavg[i]); } for(i=0;i<5;i++) { a[i]=1.0*dxy[i]/dx;//斜率 b[i]=Vavg[i]-a[i]*mavg;//截距 V0[i]=exp(b[i]); R[i]=dxy[i]/sqrt(dx*dy[i]); } for(i=0;i<5;i++) for(j=0;j<M;j++) SD[i]=SD[i]+pow((V[j][i]-a[i]*m[j]-b[i]),2); for(i=0;i<5;i++) SD[i]=sqrt(SD[i]/M); for(i=0;i<5;i++) printf("%lf\t%lf\t%lf\n",exp(b[i]),R[i],SD[i]); //水汽波段定标 //计算大气柱总光学厚度 for(i=0;i<M;i++) for(j=0;j<5;j++) V[i][j]=exp(V[i][j]); for(i=0;i<M;i++) for(j=0;j<5;j++) tao[i][j]=log(V0[j]/V[i][j])/m[i]; for(i=0;i<5;i++) { taor[i]=0.008569*pow(namta8[i],-4)*(1+0.0113*pow(namta8[i],-2)+0.00013*pow(namta8[i],-4))*1000*exp(-0.125017*0.044)/1013.25; //计算大气分子垂直光学厚度taor[i]=0.0088*pow(namta8[i],-4.05) //printf("%lf\t",taor[i]); } for(i=0;i<M;i++) for(j=0;j<5;j++) taoa8[i][j]=tao[i][j]-taor[j]; //计算气溶胶光学厚度，含水汽 //利用线性内插法计算940波段气溶胶光学厚度 for(i=0;i<M;i++) { taoa940[i]=taoa8[i][0]+(taoa8[i][2]-taoa8[i][0])*(namta8[1]-namta8[0])/(namta8[2]-namta8[0]); taoH2O[i]=taoa8[i][1]-taoa940[i]; //printf("%lf\n",taoa940[i]); } /*利用拟合法计算550nm,940nm波段气溶胶光学厚度：气溶胶光学厚度在0.4~1.1um满足Angstrom 公式：taoanamta=beta*namta.^arfa,beta为浑浊度系数,arfa反映了大气气溶胶粒子谱分布*/ //拟合得到Angstram波长指数和浑浊度系数 dx=0.0; total=0.0; for(i=0;i<M;i++) for(j=0;j<4;j++) if(j>0) taoa5[i][j]=taoa8[i][j+1]; //除936波段以外的气溶胶光学厚度 else taoa5[i][j]=taoa8[i][j]; for(i=0;i<M;i++) for(j=0;j<4;j++) taoa5[i][j]=log(taoa5[i][j]); for(i=0;i<4;i++) { if(i>0) lnnamta5[i]=log(namta8[i+1]); else lnnamta5[i]=log(namta8[i]); total=total+lnnamta5[i]; } namta5avg=total/4; for(i=0;i<M;i++) for(j=0;j<4;j++) sumtaoa5[i]=sumtaoa5[i]+taoa5[i][j]; for(i=0;i<M;i++) taoa5avg[i]=sumtaoa5[i]/4; for(i=0;i<4;i++) dx=dx+pow((lnnamta5[i]-namta5avg),2); for(i=0;i<M;i++) for(j=0;j<4;j++) dxym[i]=dxym[i]+(lnnamta5[j]-namta5avg)*(taoa5[i][j]-taoa5avg[i]); for(i=0;i<M;i++) { slope1[i]=1.0*dxym[i]/dx; //斜率 Angstram波长指数 arfa[i]=-slope1[i]; beta[i]=taoa5avg[i]-slope1[i]*namta5avg; //截距 //beta[i]=exp(beta[i]); //浑浊度系数 } for(i=0;i<M;i++) beta[i]=exp(beta[i]); //计算550nm,1050nm,940nm,520nm波段的气溶胶光学厚度 for(i=0;i<M;i++) { // taoa8[i][0]=beta[i]*pow(1.050,-arfa[i]); taoa8[i][1]=beta[i]*pow(0.936,-arfa[i]); // taoa8[i][6]=beta[i]*pow(0.520,-arfa[i]); //taoa550[i]=beta[i]*pow(0.550,-arfa[i]); } p=fopen(DEPTHPATH,"w"); if(p==NULL) { printf("can't open the file!\n"); return; } fprintf(p,"time\t1020nm\t936nm\t870nm\t670nm\t440nm\n"); for(i=0;i<M;i++) { fprintf(p,"%5.2lf\t",timal[i]); for(j=0;j<5;j++) fprintf(p,"%7.6lf\t",taoa8[i][j]); fprintf(p,"\n"); } fclose(p); //940波段分子与气溶胶光学厚度 for(i=0;i<M;i++) tao940[i]=taoa8[i][1]+taor[1]; //改进的langley法对水汽波段定标 dx=0.0; dy[1]=0.0; total=0.0; SD[1]=0.0; for(i=0;i<M;i++) { bm[i]=pow(m[i],ba); total=total+bm[i]; } bmavg=total/M; total=0.0; for(i=0;i<M;i++) { Vm[i]=log(V[i][1])+m[i]*tao940[i]; total=total+Vm[i]; } Vmavg=total/M; for(i=0;i<M;i++) { dx=dx+pow((bm[i]-bmavg),2); dy[1]=dy[1]+pow((Vm[i]-Vmavg),2); //printf("%lf\t",dx); dxy940=dxy940+(bm[i]-bmavg)*(Vm[i]-Vmavg); //printf("%lf\t",dxy940); } a[1]=1.0*dxy940/dx;//斜率 b[1]=Vmavg-a[1]*bmavg;//截距 R[1]=dxy940/sqrt(dx*dy[1]); for(i=0;i<M;i++) SD[1]=SD[1]+pow((Vm[i]-a[1]*bm[i]-b[1]),2); SD[1]=sqrt(SD[1]/M); if((p=fopen(V0RSDPATH,"a"))==NULL) { printf("can't open the file!\n"); return; } fprintf(p,"%s\n",str); fprintf(p,"V0 lnV0 R SD \n"); for(i=0;i<5;i++) fprintf(p,"%6.4lf %4.4lf %4.6lf %4.6lf\n",exp(b[i]),b[i],R[i],SD[i]); fclose(p); printf("%lf\n",exp(b[1])); if((p=fopen(ALLPATH,"w"))==NULL) { printf("打不开文件!\n"); return; } fprintf(p,"m\t1020nm\t936nm\t870nm\t670nm\t440nm\ttime\tarfa\tbeta\tbm\t936nm\n"); for(i=0;i<M;i++) { fprintf(p,"%5.4lf\t",m[i]); for(j=0;j<5;j++) fprintf(p,"%7.6lf\t",log(V[i][j])); fprintf(p,"%5.3lf\t%7.6lf\t%7.6lf\t%5.3lf\t%7.6lf\n",timal[i],arfa[i],beta[i],bm[i],Vm[i]); } fclose(p); }