/* * This file is part of tomo3d * * Copyright (C) 2002, 2003, Sebastien Judenherc <sebastien.judenherc@na.infn.it> * * 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);  }