3draytrace.c

射线追踪程序

                                } else {
                                        if (ttrefl<TT[2]) {
                                                swtch(&(param[2]),&refl,&(TT[2]),&ttrefl);
                                                sorttab(3,TT,param);
                                        } else {
                                                cont[kn]=0.5*(cent[kn]+param[2][kn]);
                                                ttcont=compT(p->nseg,X,Y,Z,T,ds,&L,cont,&(cont[NH]),NH,mod,0); L=L1;
                                                if (ttcont<TT[2]) {
                                                        swtch(&(param[2]),&cont,&(TT[2]),&ttcont);
                                                        sorttab(3,TT,param);
                                                } else {
                                                        param[1][kn]=0.5*(param[1][kn]+param[0][kn]);
                                                        TT[1]=compT(p->nseg,X,Y,Z,T,ds,&L,param[1],&(param[1][NH]),NH,mod,0); L=L1;
                                                        param[2][kn]=0.5*(param[2][kn]+param[0][kn]);
                                                        TT[2]=compT(p->nseg,X,Y,Z,T,ds,&L,param[2],&(param[2][NH]),NH,mod,0); L=L1;
                                                        sorttab(3,TT,param);
                                                }
                                        }
                                }
                                count++;
                        } while ((count<maxcount)&&(fabs(TT[2]-TT[0])>thr));
                }
                Mcnt++;
        } while ((Mcnt<30)&&(fabs(fastest-TT[0])>thr));

        compT(p->nseg, X, Y, Z, T, ds, &L, param[0], &(param[0][NH]), NH, mod, 1);
        fprintf(stderr,"-> %10.3f sec (%4d %e) L=%.2f/%.2f\n",TT[0],Mcnt,T0-TT[0],L1,L);
        fprintf(stderr,"PERIOD   NORM ( LONG,TRANS) N\n");
        fprintf(fout,"-> %10.3f sec ( %4d %e ) L= %.2f / %.2f\n",TT[0],Mcnt,T0-TT[0],L1,L);
        fprintf(fout,"PERIOD   NORM ( LONG,TRANS) N\n");
        for (i=0; i<NH; i++) {
                fprintf(stderr,"%6.2f %6.2f ( %5.2f,%5.2f ) %d %s\n",L1/(float)(i+1),
                                sqrt(param[0][i]*param[0][i]+param[0][NH+i]*param[0][NH+i]),
                                param[0][i],param[0][NH+i],i,
                                missed?"   MISSED":"");
                fprintf(fout,"%6.2f %6.2f ( %5.2f , %5.2f ) %d %s\n",L1/(float)(i+1),
                                sqrt(param[0][i]*param[0][i]+param[0][NH+i]*param[0][NH+i]),
                                param[0][i],param[0][NH+i],i,
                                missed?"   MISSED":"");
        }
        fflush(fout);

        for (i=0; i<3; i++) free(param[i]);
        free(param);
        free(cent);
        free(cont);
        free(refl);
        free(expn);
}

int main(int argc, char **argv)
{
        struct path_t * path=NULL;
        float *d,*sd;
        float *X=NULL,*Y=NULL,*Z=NULL,*T=NULL;
        FILE *f,*fr3d,*fharmo,*fgmt;
        float DS,W=20;
        int Nray,i,NAFF,NH=10;
        int NSEGMAX;
        char rayname[255],xyzname[255],r3dname[255],harname[255],gmtname[255];
        char *basename=NULL;
        struct path_t *p;
        struct mod3d_t *m3d=NULL;
        int GMT=0;
        char *refmodname=NULL;
        float P2Sconv=0.0;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit(0); }
        init_mod(NULL,0.);
        if (argc<2) goto usage;
        for (i=1; i<argc; i++) {
                if (argv[i][0]=='-') {
                        switch (argv[i][1]) {
                                case 'g' :
                                        GMT=1;
                                        break;
                                case 'n' :
                                        if (argc<=i+1) goto error;
                                        basename=argv[++i];
                                        break;
                                case 'W' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%f",&W)) goto error;
                                        W/=2.;
                                        break;
                                case 'H' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%d",&NH)) goto error;
                                        break;
                                case 'c' :
                                        if (argc<=i+2) goto error;
                                        refmodname=argv[++i];
                                        if (!sscanf(argv[++i],"%f",&P2Sconv)) goto error;
                                        break;
                                default  :
                                        goto error;
                        }
                } else goto error;
        }
        init_mod(refmodname,P2Sconv);
        if (!basename) goto usage;

        sprintf(gmtname,"%s.r3t",basename);
        sprintf(rayname,"%s.ray",basename);
        sprintf(xyzname,"%s.xyz",basename);
        sprintf(r3dname,"%s.r3d",basename);
        sprintf(harname,"%s.har",basename);

        if (NULL==(fharmo=fopen(harname,"wt"))) {
                perror(harname);
                exit(1);
        }

        if (NULL==(fgmt=fopen(gmtname,"wt"))) {
                perror(gmtname);
                exit(1);
        }

        /*m3d=readveloc(xyzname,5.0,20.0);*/
        m3d=readveloc(xyzname,20.0,W);

        fprintf(stderr,"Reading %s...",rayname);
        if ((f=fopen(rayname,"rb"))==NULL) {
                perror(rayname);
                exit(1);
        }
        fread(&Nray, sizeof(int), 1, f);
        fread(&DS, sizeof(float), 1, f);
        if ((fr3d=fopen(r3dname,"wb"))==NULL) {
                perror(r3dname);
                exit(1);
        }
        fwrite(&Nray, sizeof(int), 1, fr3d);
        fwrite(&DS, sizeof(float), 1, fr3d);
        /*path=readpath(rayname,&Nray,&DS);*/
        NAFF=Nray/10;
        if (!NAFF) NAFF=10;
        fprintf(stderr,"%d rays, DS=%f\n",Nray,DS);

        d= (float*)calloc(Nray,sizeof(float));
        sd=(float*)calloc(Nray,sizeof(float));

        /* allocation for all the ray headers (= path) */
        path=(struct path_t *)malloc(Nray*sizeof(struct path_t));
        NSEGMAX=-1;
        ptime(0,1);
        for (i=0; i<Nray; i++) {
                p=&(path[i]);
                /* read the ray header */
                fread(p,sizeof(struct path_t),1,f);
                if (p->nseg>NSEGMAX) {
                        NSEGMAX=p->nseg+50;
                        X=realloc(X,sizeof(float)*NSEGMAX);
                        Y=realloc(Y,sizeof(float)*NSEGMAX);
                        Z=realloc(Z,sizeof(float)*NSEGMAX);
                        T=realloc(T,sizeof(float)*NSEGMAX);
                        if ((!X)||(!Y)||(!Z)||(!T)) {
                                fprintf(stderr,"alloc error\n");
                                exit(1);
                        }
                }

                /* read the ray path and time along segments DS */
                fread(X,sizeof(float),p->nseg,f);
                fread(Y,sizeof(float),p->nseg,f);
                fread(Z,sizeof(float),p->nseg,f);
                fread(T,sizeof(float),p->nseg,f);

                if (GMT) {
                        int j;
                        for (j=0; j<p->nseg; j++)
                                fprintf(fgmt,"%.1f %.1f %.1f 1D\n",X[j],Y[j],Z[j]);
                        fprintf(fgmt,"> 1D\n");
                        fflush(fgmt);
                }

                findbestSIMPLEX(p,X,Y,Z,T,NH,m3d,DS,fharmo);
                fwrite(p,sizeof(struct path_t),1,fr3d);
                fwrite(X,sizeof(float),p->nseg,fr3d);
                fwrite(Y,sizeof(float),p->nseg,fr3d);
                fwrite(Z,sizeof(float),p->nseg,fr3d);
                fwrite(T,sizeof(float),p->nseg,fr3d);
                fflush(fr3d);
                if (GMT) {
                        int j;
                        for (j=0; j<p->nseg; j++)
                                fprintf(fgmt,"%.1f %.1f %.1f 3D\n",X[j],Y[j],Z[j]);
                        fprintf(fgmt,"> 3D\n");
                        fflush(fgmt);
                }
        }
        fclose(fr3d);
        fclose(fharmo);
        if (GMT) fclose(fgmt);
        ptime(0,2);
        return(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s [-H nharm(10)] [-W width(20)] [-g] [-c refmod-file P2Sconv] -n name\n",
                        argv[0]);
        return(1);
}