计算全球震相的程序

最新推荐文章于 2025-12-19 15:31:04 发布

原创最新推荐文章于 2025-12-19 15:31:04 发布 · 653 阅读

CC 4.0 BY-SA版权

文章标签：

int main(int argc,char*argv[]) {
   if(argc<3)
       {
           cout<<"Usage:"<<argv[0]<<"<phase/Psdep> <delta> [depth=0]"<<" [initrayp] [recdep] [rayperr] [deltaerr] [Nsamp]"<<endl;
           return 0;
           }

double sigma_misfit_restored = 1e-9;
int N = 50; // 需要根据实际需求设置
int max_iterations = 100;

double data=atof(argv[2])*PI/180;

double dep=R,surf=R,z4=1;

if(argc>3)dep=R-abs(atof(argv[3]));
if(argc>5)surf=R-atof(argv[5]);
if(argc>7)sigma_misfit_restored=atof(argv[7]);
if(argc>8){N=atoi(argv[8]);cout<<"num samples:"<<N<<endl;}
int ndata=1, nparams = 1,i;
int nm,code[NP],nmod;
   double z3[NP],rayp=10,rp[NMOD],dp=1;
FILE *fp=NULL;fp=fopen("iasp91.txt","r");
if(fp!=NULL)fclose(fp);else {fp=fopen("iasp91.txt","w"); for(i=0;i<NI-1;i++)fprintf(fp,"%10.3f %10.4f %10.4f %10.4f\n",mod0[i][0],mod0[i][1],mod0[i][2],mod0[i][3]); fprintf(fp,"%10.3f %10.4f %10.4f %10.4f",mod0[i][0],mod0[i][1],mod0[i][2],mod0[i][3]); fclose(fp);}
PSModel ps("iasp91.txt");code[0]=2;nm=1;
if(atof(argv[1])>0){code[0]=2;code[1]=-1;z3[0]=R-atof(argv[1]);nm=2;}
if(std::string(argv[1])==std::string("SKS")){code[0]=-1;code[1]=2;code[2]=-1;z3[0]=-2;z3[1]=-2;nm=3;}
if(std::string(argv[1])==std::string("SKiKS")){code[0]=-1;code[1]=2;code[2]=1;code[3]=2;code[4]=-1;z3[0]=-2;z3[1]=-1;z3[2]=-1;z3[3]=-2;nm=5;}
if(std::string(argv[1])==std::string("SKIKS")){code[0]=-1;code[1]=2;code[2]=2;code[3]=2;code[4]=-1;z3[0]=-2;z3[1]=-1;z3[2]=-1;z3[3]=-2;nm=5;}
if(std::string(argv[1])==std::string("PKiKP")){code[0]=1;code[1]=2;code[2]=1;code[3]=2;code[4]=1;z3[0]=-2;z3[1]=-1;z3[2]=-1;z3[3]=-2;nm=5;}
if(std::string(argv[1])==std::string("PKIKP")){code[0]=1;code[1]=2;code[2]=2;code[3]=2;code[4]=1;z3[0]=-2;z3[1]=-1;z3[2]=-1;z3[3]=-2;nm=5;}
if(std::string(argv[1])==std::string("P")){code[0]=2;nm=1;}
if(std::string(argv[1])==std::string("p")){code[0]=1;nm=1;}
if(std::string(argv[1])==std::string("S")){code[0]=-2;nm=1;}
if(std::string(argv[1])==std::string("s")){code[0]=-1;nm=1;}
if(std::string(argv[1])==std::string("Pm")){code[0]=1;code[1]=1;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("Sm")){code[0]=-1;code[1]=-1;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("Pms")){code[0]=1;code[1]=-1;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("Smp")){code[0]=-1;code[1]=1;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("Ps")){code[0]=2;code[1]=-1;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("Sp")){code[0]=-2;code[1]=1;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("PP")){code[0]=2;code[1]=2;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("SS")){code[0]=-2;code[1]=-2;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("PS")){code[0]=2;code[1]=-2;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("SP")){code[0]=-2;code[1]=2;z3[0]=-7;nm=2;}
if(std::string(argv[1])==std::string("PcP")){code[0]=1;code[1]=1;z3[0]=-2;nm=2;}
if(std::string(argv[1])==std::string("ScS")){code[0]=-1;code[1]=-1;z3[0]=-2;nm=2;}
if(std::string(argv[1])==std::string("PcS")){code[0]=1;code[1]=-1;z3[0]=-2;nm=2;}
if(std::string(argv[1])==std::string("ScP")){code[0]=-1;code[1]=1;z3[0]=-2;nm=2;}
if(std::string(argv[1])==std::string("Pn")){rayp=ps.Pn(dep,surf,-7,data);cout<<"rayp="<<rayp<<" time="<<ps.time<<" error="<<(ps.warc-data)*180/ PI<<" "<<ps.warc*180/PI<<endl;return 1;}
if(std::string(argv[1])==std::string("Sn")){rayp=ps.Sn(dep,surf,-7,data);cout<<"rayp="<<rayp<<" time="<<ps.time<<" error="<<(ps.warc-data)*180/ PI<<" "<<ps.warc*180/PI<<endl;return 2;}
if(std::string(argv[1])==std::string("Pdiff")){rayp=ps.Pn(dep,surf,2,data);cout<<"rayp="<<rayp<<" time="<<ps.time<<" error="<<(ps.warc-data)*180/ PI<<" "<<ps.warc*180/PI<<endl;return 3;}
if(std::string(argv[1])==std::string("Sdiff")){rayp=ps.Sn(dep,surf,2,data);cout<<"rayp="<<rayp<<" time="<<ps.time<<" error="<<(ps.warc-data)*180/ PI<<" "<<ps.warc*180/PI<<endl;return 3;}

VectorXd fop=params_final.colwise().mean();
for(int i=0;i<fop.size();i++)
{

ps.sp4(fop(i),dep,surf,z3,code,nm);
cout<<"rayp="<<fop(i)<<" time="<<ps.time<<" error="<<(ps.warc-d0(i))*180/ PI<<" "<<ps.warc*180/PI<<endl;

}
return 0;
}