/* * 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 * */#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-8double 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;