raytrace.c

射线追踪程序

/* trace de rai dans un modele "bloc"
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <eve.h>
#include <sismoutil.h>
#include "tomo3d.h"

clock_t __lclock=0;
#define PROGRESS(k,N) { if (clock()>__lclock+CLOCKS_PER_SEC) { fprintf(stderr,"%6d/%6d",k,N); __lclock=clock(); } }

#define RE 6371.
#define STEPALLOC 50

#define OK fprintf(stderr,"FILE %s, LINE %d\n",__FILE__,__LINE__);

double homo(double x)
{
        return(8.0);
}

int main(int argc, char **argv)
{
        struct event_t *cur=NULL;
        struct path_t *path=NULL;
        struct data_t *data=NULL;
        struct data_t *dcur=NULL;
        char *basename=NULL;
        char rayname[255],gmtname[255],layname[255],staname[255];
        char stacode[1000][5];
        struct lambert proj;
        float DS=0.5,maxdepth=-1.0;
        float VEL;
        time_t iclock;
        int i,j,MINRAY=3,Nsta=0,Nray=0,Nray0;
        int GMTdump=0,norm=0;
        FILE *f,*fray;
        char *refmodname=NULL;
        float P2Sconv=0.0;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit; }
        init_mod(NULL,0.);
        if (argc<2) goto usage;
        for (i=1; i<argc; i++) {
                if (argv[i][0]=='-') {
                        switch (argv[i][1]) {
                                case 'c' :
                                        if (argc<=i+2) goto error;
                                        refmodname=argv[++i];
                                        if (!sscanf(argv[++i],"%f",&P2Sconv)) goto error;
                                        init_mod(refmodname,P2Sconv);
                                        break;
                                case 'i' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%f",&DS)) goto error;
                                        break;
                                case 'n' :
                                        if (argc<=i+1) goto error;
                                        basename=argv[++i];
                                        break;
                                case 's' :
                                case 'h' :
                                case 'p' :
                                        if (argc<=i+1) goto error;
                                        if (data) {
                                                data->next=(struct data_t*)malloc(sizeof(struct data_t));
                                                dcur=data->next;
                                        } else dcur=data=(struct data_t*)malloc(sizeof(struct data_t));
                                        dcur->resname=argv[i+1];
                                        switch (argv[i][1]) {
                                                case 'h' : dcur->vel=homo;
                                                           fprintf(stderr,"Pdata: %s\n",dcur->resname);
                                                           break;
                                                case 'p' : dcur->vel=Pvelocity;
                                                           fprintf(stderr,"Pdata (IASP): %s\n",dcur->resname);
                                                           break;
                                                case 's' : dcur->vel=Svelocity;
                                                           fprintf(stderr,"Sdata (IASP): %s\n",dcur->resname);
                                                           break;
                                        }
                                        dcur->next=NULL;
                                        dcur->eve=NULL;
                                        i++;
                                        break;
                                case 'm' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%d",&MINRAY)) goto error;
                                        break;
                                case 'Z' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%f",&maxdepth)) goto error;
                                        break;
                                case 'g' :
                                        GMTdump=1;
                                        break;
                                case 'N' :
                                case 'e' :
                                        if (norm) goto error;
                                        norm=2;
                                        break;
                                case 'b' :
                                        if (norm) goto error;
                                        norm=1;
                                        break;
                                default  :
                                        goto error;
                        }
                } else goto error;
        }
        if (!basename) goto usage;
        if (!data) goto usage;
        if (maxdepth<0.) goto error;

        fprintf(stderr,"================= raytrace ================\n");

        sprintf(rayname,"%s.ray",basename);
        sprintf(gmtname,"%s.gmt",basename);
        sprintf(layname,"%s.lay",basename);
        sprintf(staname,"%s.sta",basename);

        f=fopen(layname,"rt");
        if (!f) {
                perror(layname); exit(0);
        }
        else {
                double lat=48.0,lon=-2.5,X,Y;
                if (fscanf(f,"%lf %lf",&lat,&lon)!=2) {
                        fprintf(stderr,"ERR: %s\n",layname);
                        exit(0);
                }
                initproj(13,lon,lat,0.0,&proj);
                geo2lam(lat,lon,&X,&Y,&proj);
                fprintf(stderr,"PROJ: %s, ",proj.name);
        }
        fprintf(stderr,"DS=%.1f, ",DS);
        fclose(f);

        fprintf(stderr,"Zmax=%.0f, ",maxdepth);
        fprintf(stderr,"MINRAY=%d\n",MINRAY);

        for (dcur=data; dcur; dcur=dcur->next) {
                fprintf(stderr,"Reading %s...",dcur->resname);
                if ((dcur->eve=readevent(dcur->resname))==NULL) exit(0);
        }
        fprintf(stderr," done\n");

        fprintf(stderr,"  Selecting events...");
        Nray=0;
        for (dcur=data; dcur; dcur=dcur->next)
                for (cur=dcur->eve; cur; cur=cur->next) {
                        if (cur->Npick<MINRAY) cur->Npick=0; else Nray+=cur->Npick;
                }
        fprintf(stderr," Cleaning/numbering events... ");
        for (dcur=data; dcur; dcur=dcur->next)
                dcur->eve=cleanevent(dcur->eve);

        if (norm>0) {
                fprintf(stderr,"Baseline normalisation... ");
                for (dcur=data; dcur; dcur=dcur->next) {
                        float locresref=0.;
                        int npck=0;
                        for (cur=dcur->eve; cur; cur=cur->next)
                                for (i=0; i<cur->Npick; i++,npck++) locresref+=cur->pick[i].res;
                        locresref/=(float)npck;
                        for (cur=dcur->eve; cur; cur=cur->next)
                                for (i=0; i<cur->Npick; i++) cur->pick[i].res-=locresref;
                        fprintf(stderr,"average residual (%s): %.3f\n",
                                        dcur->resname,locresref);
                }
        }
        if (norm>1) {
                fprintf(stderr,"Event normalisation...\n");
                for (dcur=data; dcur; dcur=dcur->next) {
                        for (cur=dcur->eve; cur; cur=cur->next) {
                                float locresref=0.;