raytraceaniso.c

射线追踪程序

  fprintf(stderr,"Selecting events...\n");
  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...\n");
  for (dcur=data; dcur; dcur=dcur->next)
      dcur->eve=cleanevent(dcur->eve);

  if (norm>0) {
     fprintf(stderr,"Baseline normalisation...\n");
     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.;
             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,"%d used\n",Nray);
  NAFF=Nray/10;
  if (!NAFF) NAFF=10;

  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);
  Nray=0;
  ptime(0,1);
  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);
                    */
                    /*inc=asin(path[Nray].p*(Pvelocity(Z)));*/
                    VEL=dcur->vel(Z);
                    _s=path[Nray].p*VEL;
                    assert(fabs(_s)<=1.0);
                    /*
                    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*_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;
  /*******************************************************************/
                    /* quel modele aniso ? */
                    /* 30km vers le sud */
                    if (path[Nray].y[ns] >
                        (-.24302981197291080152*(path[Nray].x[ns]-423.930)+2280.434)) {
                       aniso=anisoNord;
                    } else aniso=anisoSud;
                    /* NAD+CAD
                    if (path[Nray].y[ns] >
                        (-.24302981197291080152*(path[Nray].x[ns]-423.930)+2280.434)) {
                       aniso=anisoNord;
                    } else aniso=anisoSud;
                    */
                    if (path[Nray].z[ns]<150.0)
                       VEL*=perturbaniso(path[Nray].az,path[Nray].i,aniso,5.0);
  /*******************************************************************/
                    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++;
              if (! (Nray%NAFF)) fprintf(stderr,"%d ",Nray);
          } /* 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);
}