raytrace2.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
 *
 */
/* fabrique un fichier d'entree pour le dessin des rais en 3D
 * format de sortie= x y z */

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

#define RE 6370.

int main(int argc, char **argv)
{ struct block_t * blk=NULL;
  struct event_t *eve=NULL;
  struct event_t *cur=NULL;
  struct ray_t *ray=NULL;
  struct ray_t *baseray=NULL;
  struct lambert proj;
  FILE *fdes;
  int dessin=1;
  int nbb,num,lastnum,ip,Nray=0,OK;
  char stacode[1000][5];
  int Nsta=0,is,layer1=0;
  int j,i,stanum;
  int Neve=0;
  double ds=0.5;
  double dx,dy,dz,X,Y;
  /*initproj(LAMB_II_E,0.0,0.0,0.0,&proj);*/
  argv=parseoptions(&argc,argv);
  if (argc<0) { argc*=-1; exit; }
  initproj(PROJ_UTM,-2.5,0.0,0.0,&proj);
  if (argc<3) {
     fprintf(stderr,"%s <block-file> <event-file>\n",argv[0]);
     exit(0);
  }
  if ((blk=readblock(argv[1], &nbb))==NULL) exit(0);
  for (layer1=0,X=blk[0].z; layer1<nbb; layer1++) {
      if (blk[layer1].z>X+ds) break;
      X=blk[layer1].z;
  }
  /*for (i=0; i<nbb; i++) blk[i].dV*=1.0;*/
  fdes=fopen("rais.dat","wt");
  fprintf(stderr,"Layer 1 ending at %d\n",layer1);
  fprintf(stderr,"Max depth=%.1f\n",blk[nbb-1].z+blk[nbb-1].dz);
  fprintf(stderr,"Reading %s...",argv[2]);
  fflush(stderr);
  if ((eve=readevent(argv[2]))==NULL) exit(0);
  fprintf(stderr,"OK\n");
  
  fprintf(stderr,"Tracing rays ...");
  for (cur=eve; cur; cur=cur->next) {
      for (OK=0,ip=0; ip<cur->Npick; ip++) {
          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);
             fprintf(stderr,"%d %s\n",stanum,stacode[Nsta]);
             Nsta++;
          }
          geo2lam(cur->pick[ip].sta->lat,cur->pick[ip].sta->lon,&X,&Y, &proj);
          lastnum=num=blocknum(X/1000.,Y/1000.,0.0,blk,nbb,0);
          if (num<0) {       /* is the station in the target volume? */
             fprintf(stderr,"starting out!!\n");
             continue;
          }
          /*if (num<layer1) num=stanum;*/
          OK=1;
          /* yes: new ray */
          baseray=(struct ray_t *)realloc(baseray,
                              (Nray+1)*sizeof(struct ray_t));
          ray=&(baseray[Nray]);
          Nray++;
          ray->stanum=stanum;
          ray->block=num;
          /* m -> km */
          ray->x=X/1000.;
          ray->y=Y/1000.;
          ray->z=0.1;
          if (dessin) fprintf(fdes,">\n%f %f 0.0\n",ray->x,ray->y);
          /*
          fprintf(stderr,
             "Source %4.4s X=%.0f Y=%.0f: lat=%5.1f lon=%5.1f baz=%5.1f\n",
                         cur->pick[ip].sta->code,
                         ray->x,ray->y,
                         cur->pick[ip].sta->lat,
                         cur->pick[ip].sta->lon,
                         cur->pick[ip].baz);
                         */
          ray->evenum=Neve;
          ray->i=(float)asin((blk[num].V*(1./*+blk[num].dV*/)
                             *(cur->pick[ip].slow)*180.)/(M_PI*RE));
          ray->p=sin(ray->i)/(blk[num].V*(1.+blk[num].dV));
          ray->az=M_PI*(cur->pick[ip].baz)/180.;
          ray->L=0.;
          ray->res=0.0;
          /*ray->res=0;*/
          ray->next=NULL;
          while ((num=blocknum(ray->x,ray->y,ray->z,blk,nbb,lastnum))>0) {
                if (dessin) fprintf(fdes,"%f %f %f\n",ray->x,ray->y,ray->z);
                /*if (num<layer1) num=stanum;*/
                ray->i=asin(ray->p*(blk[num].V*(1./*+blk[num].dV*/)));
                dx=ds*sin(ray->az)*sin(ray->i);
                dy=ds*cos(ray->az)*sin(ray->i);
                dz=ds*cos(ray->i);
                if (dx>3*ds) fprintf(stderr,"!!! dx=%f\n",dx);
                if (dy>3*ds) fprintf(stderr,"!!! dy=%f\n",dy);
                ray->x+=dx;
                ray->y+=dy;
                ray->z+=dz;
                ray->block=num;
                if (num!=lastnum) {
                   baseray[Nray-1].res+=ray->L;
                        /*
                   fprintf(stderr,"Changement de bloc...%d (%3.1f)->%d\n",
                                  lastnum,ray->L,num);
                                  */
                   ray->next=(struct ray_t *)malloc(sizeof(struct ray_t));
                   memcpy(ray->next,ray,sizeof(struct ray_t));
                   ray->next->next=NULL;
                   ray=ray->next;
                   ray->L=(ds/(blk[num].V*(1.0+blk[num].dV))
                          -ds/blk[num].V);
                   lastnum=num;
                } else ray->L+=(ds/(blk[num].V*(1.0+blk[num].dV))
                               -ds/blk[num].V);
          }
          baseray[Nray-1].res+=ray->L;
          cur->pick[ip].res-=baseray[Nray-1].res;
          if (ray->z<blk[nbb-1].z+blk[nbb-1].dz) {
             struct ray_t *cur;
             fprintf(stderr,"Zout<Zmax Nray=%d\n",Nray-1);
             Nray--;
             while (&(baseray[Nray])) {
                   fprintf(stderr,"."); fflush(stderr);
                   for (cur=&(baseray[Nray]); cur->next; cur=cur->next) ;
                   free(cur);
                   cur=NULL;
             }

          }
          /*
          fprintf(stderr,"Out point: %f %f %f\n",ray->x,ray->y,ray->z);
          */
      }
      if (OK) Neve++;
  }
  if (dessin) fclose(fdes);
  fprintf(stderr,"OK\n");
  fprintf(stderr,"%d rays, %d stations, %d events\n",Nray, Nsta, Neve);
  writedat(stdout,eve);
  return(0);
}