block2dat.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 <sismoutil.h>
#include <math.h>
#include "tomo3d.h"

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

#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;
        double lat=48.0,lon=-2.5;
        struct lambert proj;
        char blkname[255],parname[255],numname[255],staname[255];
        char xyzname[255],geoname[255],hitname[255],layname[255];
        char gmtname[255];
        FILE *f,*flay,*fnum,*fgeo,*fxyz,*fhit,*fsta;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit(0); }
        if (argc<2) {
                fprintf(stderr,"%s name\n",argv[0]);
                exit(0);
        }
        fprintf(stderr,"================ block2dat ================\n");
        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]);
        sprintf(gmtname,"%s-blk.gmt",argv[1]);
        flay=fopen(layname,"rt");
        if (!flay) {
                perror(layname);
                exit(0);
        }
        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);

        f=fopen(gmtname,"wt");
        for (bl=0; bl<nbb; bl++) {
                fprintf(f,"> %.0f\n",blk[bl].z);
                lam2geo(1000.*(blk[bl].x),
                                1000.*(blk[bl].y),
                                &lat,&lon,&proj);
                fprintf(f,"%f %f %.0f\n",lat,lon,blk[bl].z);
                lam2geo(1000.*(blk[bl].x+blk[bl].dx),
                                1000.*(blk[bl].y),
                                &lat,&lon,&proj);
                fprintf(f,"%f %f %.0f\n",lat,lon,blk[bl].z);
                lam2geo(1000.*(blk[bl].x+blk[bl].dx),
                                1000.*(blk[bl].y+blk[bl].dy),
                                &lat,&lon,&proj);
                fprintf(f,"%f %f %.0f\n",lat,lon,blk[bl].z);
                lam2geo(1000.*(blk[bl].x),
                                1000.*(blk[bl].y+blk[bl].dy),
                                &lat,&lon,&proj);
                fprintf(f,"%f %f %.0f\n",lat,lon,blk[bl].z);
                lam2geo(1000.*(blk[bl].x),
                                1000.*(blk[bl].y),
                                &lat,&lon,&proj);
                fprintf(f,"%f %f %.0f\n",lat,lon,blk[bl].z);
        }
        fprintf(f,">");
        return(0);
}