syntheaniso.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
 *
 */
#define __tomo3d_c
#include <math.h>
#include <sismoutil.h>
#include <eigen.h>
#include <aniso.h>
#include "tomo3d.h"

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

extern float RAYAZ;
extern float RAYINC;

float ANISOTX=0.0;
float VISO=1.0;

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])/VISO)-1.0)/ANISOTX);
}

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=1.;
  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);
}

void initaniso2(float az1, float dip1, float az2, float dip2)
{ m=matrixalloc(3);
  pol=matrixalloc(3);
  prop=vectoralloc(3);
  veloc=vectoralloc(3);
  azNord=az1;
  dipNord=dip1;
  azSud=az2;
  dipSud=dip2;
  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);
  { float az,inc;
    double azd,dipd;
    float faz,fdip;
    float Vmin=1000.0,Vmax=-1000.,Vmoy=0.,V;
    int Nv=0;
    FILE *f;
    f=fopen("vaniso.bin","wb");
    fprintf(stderr,"V: min, max, iso, aniso: ");
    for (az=0.0; az<360.0; az+=1.0)
        for (inc=-1.0; inc<92.0; inc+=1.0) {
            geoph2xyz(az,inc,&(prop[0]),&(prop[1]),&(prop[2]));
            m=christoffel(aniso,prop,m);
            eigen(m,3,veloc,pol,3);
            V=sqrt(veloc[0]);
            xyz2azdip(prop[0],prop[1],prop[2],&azd,&dipd);
            faz=(float)azd;
            fdip=-(float)dipd;
            fwrite(&faz,sizeof(float),1,f);
            fwrite(&fdip,sizeof(float),1,f);
            fwrite(&V,sizeof(float),1,f);
            Vmoy+=V;
            Nv++;
            if (V<Vmin) Vmin=V;
            if (V>Vmax) Vmax=V;
        }
    Vmoy/=(float)Nv;
    fprintf(stderr,"%.3f %.3f %.3f %.2f%%\n",
            Vmin,Vmax,Vmoy,100.0*(Vmax-Vmin)/Vmoy);
    ANISOTX=100.0*(Vmax-Vmin)/Vmoy;
    VISO=Vmoy;
  }
}

void initaniso()
{ 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);
  { float az,inc;
    double azd,dipd;
    float faz,fdip;
    float Vmin=1000.0,Vmax=-1000.,Vmoy=0.,V;
    int Nv=0;
    FILE *f;
    f=fopen("vaniso.bin","wb");
    fprintf(stderr,"V: min, max, iso, aniso: ");
    for (az=0.0; az<360.0; az+=1.0)
        for (inc=-1.0; inc<92.0; inc+=1.0) {
            geoph2xyz(az,inc,&(prop[0]),&(prop[1]),&(prop[2]));
            m=christoffel(aniso,prop,m);
            eigen(m,3,veloc,pol,3);
            V=sqrt(veloc[0]);
            xyz2azdip(prop[0],prop[1],prop[2],&azd,&dipd);
            faz=(float)azd;
            fdip=-(float)dipd;
            fwrite(&faz,sizeof(float),1,f);
            fwrite(&fdip,sizeof(float),1,f);
            fwrite(&V,sizeof(float),1,f);
            Vmoy+=V;
            Nv++;
            if (V<Vmin) Vmin=V;
            if (V>Vmax) Vmax=V;
        }
    Vmoy/=(float)Nv;
    fprintf(stderr,"%.3f %.3f %.3f %.2f%%\n",
            Vmin,Vmax,Vmoy,100.0*(Vmax-Vmin)/Vmoy);
    ANISOTX=100.0*(Vmax-Vmin)/Vmoy;
    VISO=Vmoy;
  }
}

float ANISO1(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  if (z>P) return(1.0);
  A*=100.;
  /* >48N: */
  if (y<0.0) aniso=anisoNord; else return(1.0);
  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
  return(1.0);
}

float ANISO2(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  if (z>P) return(1.0);
  A*=100.;
  /* 1 bande 50km E-W dir rap plonge vers le sud */
  /* if ((y<-25.0)||(y>25.0)) return(1.0); */

  /* >48N: */
  if (y>0.0) aniso=anisoNord;
  /* <48N:*/
  else aniso=anisoSud;
  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
  return(1.0);
}

float ANISO3(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  if (z>P) return(1.0);
  A*=100.;
  /* 1 bloc 50x50km dir rap plonge vers le sud */
  if ((y<-25.0)||(y>25.0)) return(1.0);
  if ((x<-25.0)||(x>25.0)) return(1.0);

  aniso=anisoNord;

  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
  return(1.0);
}

float ANISO4(float x, float y, float z, float A, float P)
{ return(2.0-ANISO2(x,y,z,A,P)); }

float ANISO_G1(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  if (z>125.0) return(1.0);
  A*=100.;

  /* >48N: */
  if (y>0.0) {
     aniso=anisoNord;
     if (z>80.0) return(1.0);
  }
  /* <48N:*/
  else { aniso=anisoSud;
         if (z>125.0) return(1.0);
       }
  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
}

float ANISO_G2(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  if (z>125.0) return(1.0);
  A*=100.;

  /* Nord: */
  if (y>30-0.864113273563650635*x) {
     aniso=anisoNord;
     if (z>80.0) return(1.0);
  }
  /* Sud: */
  else { aniso=anisoSud;
         if (z>125.0) return(1.0);
       }
  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
}

float ANISO_G3(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  if (z>125.0) return(1.0);
  A*=100.;

  /* Nord: */
  if (y>-0.864113273563650635*x) {
     aniso=anisoNord;
     if (z>80.0) return(1.0);
  }
  /* Sud: */
  else { aniso=anisoSud;
         if (z>125.0) return(1.0);
       }
  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));
  m=christoffel(aniso,prop,m);
  eigen(m,3,veloc,pol,3);
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
}