/* * 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 <sismoutil.h>#include <math.h>#include "tomo3d.h"#define RE 6370.int main(int argc, char **argv){        struct block_t * blk=NULL;        int bl,nbb,lay=0,ipar,hit;        float pert,sd,res,lastdep;        struct lambert proj;        char blkname[255],parname[255],numname[255],staname[255];        char xyzname[255],geoname[255],hitname[255],layname[255];        FILE *f,*flay,*fnum,*fgeo,*fxyz,*fhit,*fsta;        argv=parseoptions(&argc,argv);        if (argc<0) { argc*=-1; exit; }        if (argc<2) {                fprintf(stderr,"%s name\n",argv[0]);                exit(0);        }        sprintf(hitname,"%s.hit",argv[1]);        sprintf(blkname,"%s.blk",argv[1]);        sprintf(parname,"%s.par",argv[1]);        sprintf(numname,"%s.num",argv[1]);        sprintf(xyzname,"%s.xyz",argv[1]);        sprintf(geoname,"%s.geo",argv[1]);        sprintf(layname,"%s.lay",argv[1]);        sprintf(staname,"%s.sta",argv[1]);        flay=fopen(layname,"rt");        if (!flay) { perror(layname); exit(0); }        {                double lat=48.0,lon=-2.5;                if (fscanf(flay,"%lf %lf",&lat,&lon)!=2) {                        fprintf(stderr,"ERR: %s\n",layname);                        exit(0);                }                initproj(13,lon,lat,0.0,&proj);        }        if ((blk=readblock(blkname, &nbb))==NULL) exit(0);        if ((fsta=fopen(staname,"rt"))==NULL) {                perror(staname);                exit(0);        }        if ((f=fopen(parname,"rt"))==NULL) {                perror(parname);                exit(0);        }        if ((fnum=fopen(numname,"rt"))==NULL) {                perror(numname);                exit(0);        }        if ((fxyz=fopen(xyzname,"wt"))==NULL) {                perror(xyzname);                exit(0);        }        if ((fgeo=fopen(geoname,"wt"))==NULL) {                perror(geoname);                exit(0);        }        if ((fhit=fopen(hitname,"rt"))==NULL) {                perror(hitname);                exit(0);        }        fscanf(f,"%d %f %f %f",&ipar,&pert,&sd,&res);        fscanf(fnum,"%d",&bl);        fscanf(fhit,"%d",&hit);        lastdep=-10.;        lay=-1;        while (!feof(f)) {                double lat,lon;                if (lastdep<blk[bl].z-0.01) {                        lay++;                        fprintf(stderr,"layer %d (%d,%d)\n",lay,bl,ipar);                }                if (lay>0) {                        fprintf(fxyz,"%10f %10f %10f %10f %10f layer%d %d %10f\n",                                        blk[bl].x+0.5*blk[bl].dx,                                        blk[bl].y+0.5*blk[bl].dy,                                        blk[bl].z+0.5*blk[bl].dz,                                        -pert*100.0,fabs(sd*100.),lay,hit,res);                        lam2geo(1000.*(blk[bl].x+0.5*blk[bl].dx),                                        1000.*(blk[bl].y+0.5*blk[bl].dy),                                        &lat,&lon,&proj);                        fprintf(fgeo,"%10e %10e %10e %10e %10e layer%d %d %10f\n",                                        lon,lat,                                        blk[bl].z+0.5*blk[bl].dz,                                        -pert*100.0,fabs(sd*100.0),lay,hit,res);                } else {                        char sta[255];                        double X,Y;                        fscanf(fsta,"%*d %s %lf %lf",sta,&lat,&lon);                        geo2lam(lat,lon,&X,&Y,&proj);                        X/=1000.;                        Y/=1000.;                        fprintf(fxyz,"%10f %10f %10f %10f %10f layer%d %d %10f\n",                                        X,Y,                                        blk[bl].z+0.5*blk[bl].dz,                                        pert*100.0,sd,lay,hit,res);                        fprintf(fgeo,"%10e %10e %10e %10e %10e layer%d %d %10f\n",                                        lon,lat,                                        blk[bl].z+0.5*blk[bl].dz,                                        -pert*100.0,fabs(sd*100.0),lay,hit,res);                }                lastdep=blk[bl].z;                fscanf(f,"%d %f %f %f",&ipar,&pert,&sd,&res);                fscanf(fnum,"%d",&bl);                fscanf(fhit,"%d",&hit);        }        fclose(f);        fclose(fnum);        fclose(fxyz);        fclose(fgeo);        fclose(fhit);        return(0);}