tomo3d.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
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <time.h>
#include <sismoutil.h>
#include <eve.h>
#include <math.h>

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

#define RE 6371.0

#include "tomo3d.h"

clock_t CC[10];
int Ncc=0;
int inter_cache=0;
int guess_blocknum_shift=0;

struct model_t *refmod=NULL;

#define EPS 1.0e-8

double PvelIASP_corr(double Z)
{
        return(PvelIASP(Z)*interp(refmod,RE-Z));
}

double SvelIASP_corr(double Z)
{
        float c;
        c=interp(refmod,RE-Z);
        return(SvelIASP(Z)*(1.0+2.0*(c-1.0)));
}

void init_mod(char *fn,float normv)
{
        if (fn) {
                refmod=readmodel(fn,normv);
        }
        if (!refmod) {
                Pvelocity=PvelIASP;
                Svelocity=SvelIASP;
                fprintf(stderr,"Using IASP only\n");
        } else {
                Pvelocity=PvelIASP_corr;
                Svelocity=SvelIASP_corr;
                fprintf(stderr,"Using IASP locally corrected\n");
        }
}

struct model_t * readmodel(char *fn,float normv)
{
        FILE *f;
        float r,v;
        struct model_t *tmp=NULL;
        f=fopen(fn,"rt");
        if (!f) { perror(fn); return(NULL); }
        tmp=(struct model_t*)calloc(1,sizeof(struct model_t));
        while (fscanf(f,"%f %f",&r,&v)==2) {
                tmp->r=(float*)realloc(tmp->r,(tmp->N+1)*sizeof(float));
                tmp->v=(float*)realloc(tmp->v,(tmp->N+1)*sizeof(float));
                tmp->r[tmp->N]=r;
                tmp->v[tmp->N]=v*normv;
                fprintf(stderr,"%d: %f %f\n",tmp->N,tmp->r[tmp->N],tmp->v[tmp->N]);
                tmp->N++;
        }
        fclose(f);
        fprintf(stderr,"%s: %d\n",fn,tmp->N);
        return(tmp);
}

float interp(struct model_t* mod,float r)
{ int i;
  float v,v1,v2,r1,r2;
  if (!mod) return(1.0);
  for (i=inter_cache; i<mod->N-1; i++) {
      if (i<0) i=0;
      r1=mod->r[i];
      r2=mod->r[i+1];
      if (((r1-0.005<=r)&&(r2+0.005>=r))||((r1+0.005>=r)&&(r2-0.005<=r))) {
         v1=mod->v[i];
         v2=mod->v[i+1];
         if (fabs(r2-r1)<0.01) return(0.5*(v1+v2));
         v=v1+(r-r1)*(v2-v1)/(r2-r1);
         /*fprintf(stderr,"OK %d %d\n",inter_cache,i);*/
         inter_cache=i;
         return(1.0+v);
      }
  }
  /*fprintf(stderr,"R1 %d %d\n",inter_cache,i);*/
  for (i=0; i<mod->N-1; i++) {
      r1=mod->r[i];
      r2=mod->r[i+1];
      if (((r1<=r)&&(r2>=r))||((r1>=r)&&(r2<=r))) {
         v1=mod->v[i];
         v2=mod->v[i+1];
         v=v1+(r-r1)*(v2-v1)/(r2-r1);
         inter_cache=i;
         /*fprintf(stderr,"%d %d\n",inter_cache,i);*/
         return(1.0+v);
      }
  }
  if (r<6300.0) fprintf(stderr,"interp error, r=%f\n",r);
  return(1.0);
}

int blkcmp(const void *p1, const void *p2)
{ struct block_t * b1=(struct block_t *)p1;
  struct block_t * b2=(struct block_t *)p2;
  int r=0;
  if (b1->z-EPS<b2->z) r-=8; else if (b1->z+EPS>b2->z) r+=8;
  if (b1->x-EPS<b2->x) r-=4; else if (b1->x+EPS>b2->x) r+=4;
  if (b1->y-EPS<b2->y) r-=2; else if (b1->y+EPS>b2->y) r+=2;
  return(r);
}

void ptime(int ncc, int comm)
{ switch (comm) {
         case -1: for (ncc=0; ncc<10; ncc++) CC[ncc]=0;
                  break;
         case  1: CC[ncc]=clock();
                  break;
         case  2: fprintf(stderr,"[%.2f sec]          ",
                          (float)(clock()-CC[ncc])/(float)CLOCKS_PER_SEC);
                  break;
  }
} 

void printlayers(FILE *f, struct layer_t *lay, int nbl)
{ int i;
  for (i=0; i<nbl; i++)
      fprintf(f,"%5.1f %5.1f  %4d %5.1f %4d %5.1f  %5.1f  %5.2f  %.2e\n",
                lay[i].x, lay[i].y, lay[i].Nx, lay[i].dx,
                lay[i].Ny, lay[i].dy, lay[i].dz, lay[i].V,lay[i].damp);
}

struct block_t *readblock(char *file, int *nbb)
{ FILE *f=fopen(file,"rt");
  struct block_t *tmp=NULL;
  struct block_t bid;
  int nl=0;
  float lz=0.;
  guess_blocknum_shift=0; /* nb moyen de block par couche */
  if (!f) {
     perror(file);
     return(NULL);
  }
  *nbb=0;
  fprintf(stderr,"Reading %s...",file);
  fflush(stderr);
  fscanf(f,"%lf %lf %lf %f %f %f %f %f %f",
         &bid.x, &bid.y,&bid.z,&bid.dx,&bid.dy,&bid.dz,
         &bid.V,&bid.dV,&bid.damp);
  lz=bid.z;
  while (!feof(f)) {
        tmp=(struct block_t*)realloc(tmp,(1+*nbb)*sizeof(struct block_t));
        memcpy(&(tmp[*nbb]),&bid,sizeof(struct block_t));
        (*nbb)++;
        fscanf(f,"%lf %lf %lf %f %f %f %f %f %f",
               &bid.x, &bid.y,&bid.z,&bid.dx,&bid.dy,&bid.dz,
               &bid.V,&bid.dV,&bid.damp);
        if (bid.z>lz+0.1) {
                nl++;
                lz=bid.z;
        }
  }
  nl++;
  guess_blocknum_shift=(*nbb)/nl;
  fclose(f);
  fprintf(stderr,"%d (Nlay=%d, avg=%d)\n",*nbb,nl,guess_blocknum_shift);
  return(tmp);
}

struct block_t *readlayers(char *file, int *nbb,
                           struct layer_t **layer, int *nbl)
{ FILE *f;
  double X, Y, Z=0.;
  float dx,dy,dz,damp,lat,lon;
  int Nx,Ny,i,j,k=0;
  float V;
  struct block_t *tmp=NULL;
  struct layer_t *lay=NULL;
  struct lambert proj;
  if ((f=fopen(file,"rt"))==NULL) {
     perror(file);
     return(NULL);
  }
  *nbl=0;
  if (fscanf(f,"%f %f",&lat,&lon)!=2) return(NULL);
  initproj(13,lon,lat,0.0,&proj);
  geo2lam(lat,lon,&X,&Y,&proj);
  X/=1000.;
  Y/=1000.;
  fprintf(stderr,"PROJ: %d lon=%.3f X=%.0f Y=%.0f\n",13,lon,X,Y);
  fscanf(f,"%d %f %d %f %f %f %f",
         &Nx,&dx,&Ny,&dy,&dz,&V,&damp);
  while (!feof(f)) {
        double minx, miny;
        minx=X-dx*0.5*(double)Nx;
        miny=Y-dy*0.5*(double)Ny;
        lay=(struct layer_t *)realloc(lay,((*nbl)+1)*sizeof(struct layer_t));
        lay[*nbl].x=X;   lay[*nbl].y=Y;   lay[*nbl].z=Z;
        lay[*nbl].dx=dx; lay[*nbl].dy=dy; lay[*nbl].dz=dz;
        lay[*nbl].Nx=Nx; lay[*nbl].Ny=Ny; lay[*nbl].V=V;
        lay[*nbl].damp=damp;

        tmp=(struct block_t *)realloc(tmp,(k+Nx*Ny)*sizeof(struct block_t));
        for (i=0; i<Nx; i++)
            for (j=0; j<Ny; j++) {
                tmp[k].x=minx+dx*(double)i;
                tmp[k].y=miny+dy*(double)j;
                tmp[k].z=Z;
                tmp[k].dx=dx;
                tmp[k].dy=dy;
                tmp[k].dz=dz;
                tmp[k].V=V;
                tmp[k].damp=damp;
                tmp[k].dV=0.0;
                k++;
            }
        Z+=tmp[k-1].dz;
        (*nbl)++;
        fscanf(f,"%d %f %d %f %f %f %f",
                  &Nx,&dx,&Ny,&dy,&dz,&V,&damp);
  }
  *nbb=k;
  *layer=lay;
  fclose(f);
  return(tmp);
}

void printblocks(FILE *f, struct block_t *blk, int nbl)
{ int i;

  for (i=0; i<nbl; i++)
      fprintf(f,"%15.5f %15.5f %15.5f   %15.8e %15.8e %15.8e    %5.2f  %8.4f %5e\n",
              blk[i].x,blk[i].y,blk[i].z,blk[i].dx,blk[i].dy,blk[i].dz,
              blk[i].V,blk[i].dV,blk[i].damp);
}

struct path_t * readpath(char *fn,int *N,float *DS)
{ FILE *f=fopen(fn,"rb");
  struct path_t *path=NULL;
  int i;
  if (!f) {
     perror(fn);
     return(NULL);
  }
  fread(N, sizeof(int), 1, f);
  fread(DS, sizeof(float), 1, f);
  path=(struct path_t *)malloc(sizeof(struct path_t)*(*N));
  for (i=0; i<*N; i++) {
      fread(&(path[i]),sizeof(struct path_t),1,f);
      path[i].x=malloc(sizeof(float)*path[i].nseg);
      path[i].y=malloc(sizeof(float)*path[i].nseg);
      path[i].z=malloc(sizeof(float)*path[i].nseg);
      fread(path[i].x,sizeof(float),path[i].nseg,f);
      fread(path[i].y,sizeof(float),path[i].nseg,f);
      fread(path[i].z,sizeof(float),path[i].nseg,f);
  }
  fclose(f);
  return(path);
}

int writepath(char *fn,struct path_t *path, int N,float DS)
{ FILE *f=fopen(fn,"wb");
  char *buf=NULL;
  int i;