synthetrace.c

射线追踪程序

/*
 * This file is part of tomo3d
 *
 * Copyright (C) 2002, 2003, Sebastien Judenherc <[e-mail]>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 * USA
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <eve.h>
#include <sismoutil.h>
#include "tomo3d.h"
#include "synthe.h"

#define RE 6371.
#define STEPALLOC 50

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

float RAYAZ,RAYINC;

double Grand(sig,M)
double sig;
double M;
{ double F1,F2,R;
  double NORM;
  if (sig<0.0) return(0.);
  NORM=(double)RAND_MAX;
  do { R=(double)rand();
       F1=1.0*(2.0*R/NORM-1.0);
       R=(double)rand();
       F2=1.0*(2.0*R/NORM-1.0);
       R=F1*F1+F2*F2;
  } while ((R<1.0e-5)||(R>1.0));
  R=sqrt(-2.0*log(R)/R);
  return(M+sig*F2*R);
}

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;
  struct lambert proj;
  float DS=0.5,maxdepth=-1.0,dinc;
  float VEL;
  float signoise=-1.0;
  int i,Nray=0,NAFF;
  int KEEPRES=0;
  char layname[255];
  char *basename=NULL;
  FILE *f;
  struct synthe_t *Vsynthe=NULL;
  int MODEL=0;
  char *refmodname=NULL;
  float P2Sconv;
  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 'c' : if (argc<=i+2) goto error;
                            refmodname=argv[++i];
                            if (!sscanf(argv[++i],"%f",&P2Sconv)) goto error;
                            break;
                            init_mod(refmodname,P2Sconv);
                 case 'm' : MODEL=1;
                            if (argc<=i+1) goto error;
                            if (!sscanf(argv[++i],"%f",&dinc)) goto error;
                            break;
                 case 'F' : if (argc<=i+3) goto error;
                            Vsynthe=parseVfunc(argv[i+1],argv[i+2],argv[i+3]);
                            if (!Vsynthe) goto error;
                            i+=3;
                            break;
                 case 'n' : if (argc<=i+1) goto error;
                            basename=argv[++i];
                            break;
                 case 'k' : KEEPRES=1;
                            break;
                 case 'i' : if (argc<=i+1) goto error;
                            if (!sscanf(argv[++i],"%f",&DS)) goto error;
                            break;
                 case 's' :
                 case 'p' : if (data) goto error;
                            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 'p' : dcur->vel=PvelIASP;
                                              fprintf(stderr,"Pdata: %s\n",dcur->resname);
                                              break;
                                   case 's' : dcur->vel=SvelIASP;
                                              fprintf(stderr,"Sdata: %s\n",dcur->resname);
                                              break;
                            }
                            dcur->next=NULL;
                            dcur->eve=NULL;
                            i++;
                            break;
                 case 'Z' : if (argc<=i+1) goto error;
                            if (!sscanf(argv[++i],"%f",&maxdepth)) goto error;
                            break;
                 case 'g' : if (argc<=i+1) goto error;
                            if (!sscanf(argv[++i],"%f",&signoise)) goto error;
                            fprintf(stderr,"Adding noise, sigma=%f\n",signoise);
                            break;
                 default  : goto error;
         }
      } else goto error;
  }
  if (!basename) goto usage;
  if (maxdepth<0.) goto error;
  sprintf(layname,"%s.lay",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\n",proj.name);
  }

  if (MODEL) {
     float lat,lon,Z,dV,ER,res;
     double X,Y;
     int N;
     char geoname[255];
     char lname[50];

     sprintf(geoname,"%s.geo",basename);
     if (NULL==(f=fopen(geoname,"rt"))) {
        perror(geoname);
        exit(1);
     }
     while (!feof(f)) {
           fscanf(f,"%f %f %f %f %f %s %d %f",
                    &lon,&lat,&Z,&dV,&ER,lname,&N,&res);

           geo2lam(lat,lon,&X,&Y,&proj);
           X/=1000.0;
           Y/=1000.0;

           dV=100.0*(1.0-compVsyn(Vsynthe,X,Y,Z));

           printf("%f %f %f %f %f %s %d %f\n",
                  lon,lat,Z,dV,ER,lname,N,res);
     }
     exit(0);
  }

  if (!data) goto usage;
  fprintf(stderr,"DS=%f\n",DS);
  fclose(f);

  fprintf(stderr,"Zmax=%f\n",maxdepth);

  Nray=0;
  for (dcur=data; dcur; dcur=dcur->next) {
      fprintf(stderr,"Reading %s...",dcur->resname);
      if ((dcur->eve=readevent(dcur->resname))==NULL) exit(0);
      for (cur=dcur->eve; cur; cur=cur->next) Nray+=cur->Npick;
  }
  fprintf(stderr,"%d used\n",Nray);
  NAFF=Nray/10;
  if (!NAFF) NAFF=10;

  fprintf(stderr,"Tracing rays ...");
  Nray=0;
  ptime(0,1);
  for (dcur=data; dcur; dcur=dcur->next) {
      for (cur=dcur->eve; cur; cur=cur->next) {
          int ip,ns;
          int Nalloc;
          double dx,dy,dz,Z;
          double _s;
          for (ip=0; ip<cur->Npick; ip++) {
              double X,Y;
              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]);
              RAYAZ=path[Nray].az=cur->pick[ip].baz*M_PI/180.;
              path[Nray].p=cur->pick[ip].slow*180./(M_PI*RE);
              RAYINC=path[Nray].i=asin(VEL*path[Nray].p);
              path[Nray].res=cur->pick[ip].res;
              path[Nray].nseg=1;
              path[Nray].T[0]=DS/(compVsyn(Vsynthe,path[Nray].x[0],
                                                   path[Nray].y[0],
                                                   path[Nray].z[0])*VEL)-DS/VEL;
              while (Z<maxdepth) {
                    VEL=dcur->vel(Z);
                    _s=path[Nray].p*VEL;
                    if (fabs(_s)>1.0) {
                       fprintf(stderr,"sin(i)=%f!!!\n",_s);
                       fprintf(stderr,"Z=%f p=%f V=%f\n",
                               Z,path[Nray].p, VEL);
                       exit(0);
                    }
                    dx=DS*sin(path[Nray].az)*_s;
                    dy=DS*cos(path[Nray].az)*_s;
                    dz=DS*sqrt(1.-_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/(compVsyn(Vsynthe,path[Nray].x[ns],
                                                          path[Nray].y[ns],
                                                          path[Nray].z[ns])*VEL)-DS/VEL;
                    Z=path[Nray].z[ns];
                    path[Nray].nseg++;
              } /* while depth... */

              if (!KEEPRES) cur->pick[ip].res=0.0;

              for (ns=0; ns<path[Nray].nseg; ns++)
                  if (KEEPRES) cur->pick[ip].res-=path[Nray].T[ns];
                     else cur->pick[ip].res+=path[Nray].T[ns];

              cur->pick[ip].res+=Grand(signoise,0.0);
              cur->pick[ip].obs=cur->pick[ip].cal+cur->pick[ip].res;
              Nray++;
              if (! (Nray%NAFF)) fprintf(stderr,"%d ",Nray);
          } /* for ip... */
      } /* for cur... */
  }
  ptime(0,2);
  fprintf(stderr,"\n");
  writedat(stdout,data->eve);
  exit(0);
error:
  fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
  fprintf(stderr,"%s -n name [-k] [-g sigma_noise] [-c refmod-file P2Sconv] [-m] -Z maxdepth [-p|s resfile] [-i inc] -F funcname P1 P2\n",argv[0]);
  fprintf(stderr,"Available functions:\n");
  printfsynlist();
  return(0);
}