raytraceaniso.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
 *
 */
/* trace de rai dans un modele "bloc" et eciture du fichier
 *  * system-name (binaire)
 *   * */

/* revblk[numero de bloc vrai]=indice de G */
/* numblk[indice de G]=numero de bloc vrai */

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

#define RE 6371.
#define STEPALLOC 50

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

scalar alphaSud, betaSud, gammaSud;
scalar alphaNord, betaNord, gammaNord;
scalar azNord=230, dipNord=30;
scalar azSud=20., dipSud=0.;
matrix ANord=NULL;
matrix ASud=NULL;
matrix m;
matrix pol;
vector prop;
vector veloc;
struct param_t *param=NULL;
struct aniso_t *aniso=NULL;
struct aniso_t *anisoNord=NULL;
struct aniso_t *anisoSud=NULL;

float perturbaniso(float az,float i,struct aniso_t *aniso,float pc)
{ geoph2xyz(az*180.0/M_PI,i*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+pc*((sqrt(veloc[0])/8.21)-1.0)/5.0);
}

scalar euler(az0,dip0,alpha, beta, gamma)
scalar az0,dip0,*alpha,*beta,*gamma;
{ scalar am=0.,bm=0.,gm=0.,a=0.,b=0.,g=0.,MF=1000000.0,mf=1000000.0;
  matrix A=NULL;
  float PAS=3.;
  vector v=vectoralloc(3);
  vector v1=vectoralloc(3);
  vector v0=vectoralloc(3);
  v[0]=0.;
  v[1]=1.;
  v[2]=0.;
  geoph2xyz(az0,90.-dip0,&(v0[0]),&(v0[1]),&(v0[2]));
  for (a=0.; a<180.; a+=PAS)
      for (b=0.; b<90.; b+=PAS)
          for (g=0.; g<90.; g+=PAS) {
              A=arot(a,b,g,A);
              rotV(v,A,v1);
              mf=acos(fabs((v0[0]*v1[0])+(v0[1]*v1[1])+(v0[2]*v1[2])))
                 *180./M_PI;
              if (mf<MF) { MF=mf; am=a; bm=b; gm=g; }
          }
  A=arot(am,bm,gm,A);
  rotV(v,A,v1);
  v[0]=am; v[1]=bm; v[2]=gm;
  printvector(stderr, "%f ", v0, 3);
  printvector(stderr, "%f ", v1, 3);
  printvector(stderr, "%f ", v, 3);
  *alpha=am;
  *beta=bm;
  *gamma=gm;
  free(v1);
  free(v0);
  free(v);
  matrixfree(A,3);
  return(MF);
}

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;
  int i,j,MINRAY=3,Nsta=0,Nray=0,NAFF;
  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(0); }
  init_mod(NULL,0.);

  /************** parametres anisotropie ****************************/
  m=matrixalloc(3);
  pol=matrixalloc(3);
  prop=vectoralloc(3);
  veloc=vectoralloc(3);
  fprintf(stderr,"MFrot=%f\n",
                euler(azNord,dipNord,&alphaNord, &betaNord, &gammaNord));
  fprintf(stderr,"MFrot=%f\n",
                euler(azSud,dipSud,&alphaSud, &betaSud, &gammaSud));
  if ((aniso=read_stiffness_hexa("hexa.dat"))==NULL) exit(0);
  anisoNord=(struct aniso_t *)malloc(sizeof(struct aniso_t));
  memcpy(anisoNord,aniso,sizeof(struct aniso_t));
  anisoSud=(struct aniso_t *)malloc(sizeof(struct aniso_t));
  memcpy(anisoSud,aniso,sizeof(struct aniso_t));
  ANord=arot(alphaNord,betaNord,gammaNord,ANord);
  ASud=arot(alphaSud,betaSud,gammaSud,ASud);
  rotC(aniso->C,ANord,anisoNord->C);
  rotC(aniso->C,ASud,anisoSud->C);
  /*******************************************************************/

  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=Pvelocity;
                                              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;

  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\n",proj.name);
  }
  fprintf(stderr,"DS=%f\n",DS);
  fclose(f);

  fprintf(stderr,"Zmax=%f\n",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);
  }