raytrace.c

射线追踪程序

                                for (i=0; i<cur->Npick; i++) locresref+=cur->pick[i].res;
                                locresref/=(float)cur->Npick;
                                for (i=0; i<cur->Npick; i++)
                                        cur->pick[i].res-=locresref;
                        }
                }
        }

        fprintf(stderr,"RAYS: %d used\n",Nray);

        fray=fopen(rayname,"wb");
        if (!fray) {
                perror(NULL);
                exit(1);
        }
        fwrite(&Nray, sizeof(int), 1, fray);
        fwrite(&DS, sizeof(float), 1, fray);
        fprintf(stderr,"Tracing rays (%s) ...",rayname);
        Nray0=Nray;
        Nray=0;
        ptime(0,1);
        iclock=time(NULL);
        for (dcur=data; dcur; dcur=dcur->next) {
                /* une serie de stations par type de residus */
                for (i=0; i<1000; i++) memset(stacode[i],0,5);
                for (cur=dcur->eve; cur; cur=cur->next) {
                        int ip,is,stanum,ns;
                        int Nalloc;
                        double dx,dy,dz,Z;
                        double _s;
                        if (cur->Npick<MINRAY) continue;
                        for (ip=0; ip<cur->Npick; ip++) {
                                double X,Y;
                                stanum=-1;
                                for (is=0; is<Nsta; is++) {
                                        if (!strcmp(cur->pick[ip].sta->code,stacode[is])) {
                                                stanum=is;
                                                break;
                                        }
                                }
                                if (stanum<0) {
                                        stanum=Nsta;
                                        strcpy(stacode[Nsta],cur->pick[ip].sta->code);
                                        Nsta++;
                                }
                                path=(struct path_t*)realloc(path,(Nray+1)*sizeof(struct path_t));
                                if (!path) {
                                        fprintf(stderr,"alloc error\n");
                                        exit(1);
                                }
                                Nalloc=STEPALLOC;
                                path[Nray].x=(float*)malloc(Nalloc*sizeof(float));
                                path[Nray].y=(float*)malloc(Nalloc*sizeof(float));
                                path[Nray].z=(float*)malloc(Nalloc*sizeof(float));
                                path[Nray].T=(float*)malloc(Nalloc*sizeof(float));
                                if ((!path[Nray].x)||(!path[Nray].y)||
                                                (!path[Nray].z)||(!path[Nray].T)) {
                                        fprintf(stderr,"alloc error\n");
                                        exit(1);
                                }
                                geo2lam(cur->pick[ip].sta->lat,
                                                cur->pick[ip].sta->lon,&X,&Y,&proj);
                                path[Nray].x[0]=X/1000.;
                                path[Nray].y[0]=Y/1000.;
                                path[Nray].z[0]=Z=0.0001;

                                VEL=dcur->vel(path[Nray].z[0]);
                                path[Nray].stanum=stanum;
                                path[Nray].evenum=cur->eve->nid;
                                path[Nray].az=cur->pick[ip].baz*M_PI/180.;
                                path[Nray].p=cur->pick[ip].slow*180./(M_PI*RE);
                                assert(fabs(VEL*path[Nray].p)<1.0);
                                path[Nray].i=asin(VEL*path[Nray].p);
                                path[Nray].res=cur->pick[ip].res;
                                path[Nray].nseg=1;
                                path[Nray].T[0]=DS/VEL;
                                /*
                                 * fprintf(stderr,"x0=%f y0=%f i=%f az=%f\n",
                                 * path[Nray].x[0],
                                 * path[Nray].y[0],
                                 * path[Nray].i*180./M_PI,
                                 * path[Nray].az*180./M_PI);
                                 * */
                                while (Z<maxdepth) {
                                        /*
                                         * fprintf(stderr,"%f %f %f/%f\n",
                                         * path[Nray].x[path[Nray].nseg-1],
                                         * path[Nray].y[path[Nray].nseg-1],
                                         * path[Nray].z[path[Nray].nseg-1],
                                         * Z,maxdepth);
                                         * */
                                        VEL=dcur->vel(Z);
                                        _s=path[Nray].p*VEL;
                                        assert(fabs(_s)<=1.0);

                                        dx=DS*sin(path[Nray].az)*_s;
                                        dy=DS*cos(path[Nray].az)*_s;
                                        dz=DS*sqrt(1.-_s*_s);
                                        if (dx>3.0*DS) fprintf(stderr,"!!! dx=%f\n",dx);
                                        if (dy>3.0*DS) fprintf(stderr,"!!! dy=%f\n",dy);
                                        if (path[Nray].nseg>Nalloc-5) {
                                                Nalloc+=STEPALLOC;
                                                path[Nray].x=
                                                        (float*)realloc(path[Nray].x,sizeof(float)*Nalloc);
                                                path[Nray].y=
                                                        (float*)realloc(path[Nray].y,sizeof(float)*Nalloc);
                                                path[Nray].z=
                                                        (float*)realloc(path[Nray].z,sizeof(float)*Nalloc);
                                                path[Nray].T=
                                                        (float*)realloc(path[Nray].T,sizeof(float)*Nalloc);
                                                if ((!path[Nray].x)||(!path[Nray].y)||
                                                                (!path[Nray].z)||(!path[Nray].T)) {
                                                        fprintf(stderr,"alloc error\n");
                                                        exit(1);
                                                }
                                        }
                                        ns=path[Nray].nseg;
                                        path[Nray].x[ns]=path[Nray].x[ns-1]+dx;
                                        path[Nray].y[ns]=path[Nray].y[ns-1]+dy;
                                        path[Nray].z[ns]=path[Nray].z[ns-1]+dz;
                                        path[Nray].T[ns]=DS/VEL;
                                        Z=path[Nray].z[ns];
                                        path[Nray].nseg++;
                                } /* while depth... */
                                fwrite(&(path[Nray]),sizeof(struct path_t),1,fray);
                                fwrite(path[Nray].x,sizeof(float),path[Nray].nseg,fray);
                                fwrite(path[Nray].y,sizeof(float),path[Nray].nseg,fray);
                                fwrite(path[Nray].z,sizeof(float),path[Nray].nseg,fray);
                                fwrite(path[Nray].T,sizeof(float),path[Nray].nseg,fray);
                                Nray++;
                                PROGRESS(Nray0-Nray,Nray0);
                        } /* for ip... */
                } /* for cur... */
                fclose(fray);
        }
        ptime(0,2);
        f=fopen(staname,"wt");
        if (!f) {
                perror(staname);
                exit(0);
        }
        for (i=0; i<Nsta; i++) {
                for (dcur=data; dcur; dcur=dcur->next) {
                        for (cur=dcur->eve; cur; cur=cur->next) {
                                int ip;
                                for (ip=0; ip<cur->Npick; ip++)
                                        if (!strcmp(cur->pick[ip].sta->code,stacode[i])) {
                                                fprintf(f,"%d %s %f %f\n",i,
                                                                cur->pick[ip].sta->code,
                                                                cur->pick[ip].sta->lat,
                                                                cur->pick[ip].sta->lon);
                                                i++;
                                        }
                        }
                }
        }
        fclose(f);
        /*writepath(rayname,path,Nray,DS);*/
        fprintf(stderr,"\n");
        if (!GMTdump) exit(0);
        fprintf(stderr,"Writing %s\n",gmtname);
        f=fopen(gmtname,"wt");
        for (i=0; i<Nray; i++) {
                fprintf(f,">\n");
                for (j=0; j<path[i].nseg; j++)
                        fprintf(f,"%f %f %f\n",path[i].x[j],path[i].y[j],path[i].z[j]);
        }
        fprintf(f,">\n");
        fclose(f);
        exit(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s [-g] -Z maxdepth -n name [-p Presfile] [-p ...] [-s Sresfile] [-s ...] [-b|e] [-i inc] [-m minray]\n",argv[0]);
        return(0);
}