crossgeo.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 <stdio.h>
#include <stdlib.h>
#include "tomo3d.h"

#define RE 6370.

struct pt_t {
        float X,Y,Z;
        float pert,res,sd,hit;
        int lay;
};

int main(int argc, char **argv)
{
        struct pt_t *pt=NULL;
        int lay=0,i,Nblk=0,j,Nit=0;
        double pert,sd,res,sum,dhit;
        char slay[255];
        float z,Z,DX,DY,L=0.;
        float minlon=-1.0,maxlon=-1.0,minlat=-1.0,maxlat=-1.0;
        double lat,lon,X,Y,minX,maxX,minY,maxY,H,depth;
        float wfact=-1.0,width=-1.0,hinc=-1.0,vinc=-1.0;
        struct lambert proj;
        char *basename=NULL;
        char geoname[255],layname[255];
        FILE *flay,*fgeo;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit(0); }
        for (i=1; i<argc; i++) {
                if (argv[i][0]=='-') {
                        switch (argv[i][1]) {
                                case 'n' :
                                        if (argc<=i+1) goto error;
                                        basename=argv[++i];
                                        break;
                                case 'i' :
                                        if (i+1>=argc) goto error;
                                        if (!sscanf(argv[++i],"%f",&hinc)) goto error;
                                        break;
                                case 'z' :
                                        if (i+1>=argc) goto error;
                                        if (!sscanf(argv[++i],"%f",&vinc)) goto error;
                                        break;
                                case 'f' :
                                        if (i+4>=argc) goto error;
                                        Nit=1;
                                        if (!sscanf(argv[++i],"%f",&minlat)) goto error;
                                        if (!sscanf(argv[++i],"%f",&minlon)) goto error;
                                        if (!sscanf(argv[++i],"%f",&maxlat)) goto error;
                                        if (!sscanf(argv[++i],"%f",&maxlon)) goto error;
                                        break;
                                case 'w' :
                                        if (i+1>=argc) goto error;
                                        if (!sscanf(argv[++i],"%f",&wfact)) goto error;
                                        break;
                                default:   goto usage;
                        }
                }
        }
        if (Nit<0) goto usage;
        if (!basename) goto usage;
        if (hinc<0.) goto usage;
        if (vinc<0.) goto usage;
        if (wfact<0.) goto usage;

        sprintf(geoname,"%s.geo",basename);
        sprintf(layname,"%s.lay",basename);
        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 ((fgeo=fopen(geoname,"rt"))==NULL) {
                perror(geoname);
                exit(0);
        }

        geo2lam(minlat,minlon,&minX,&minY,&proj);
        geo2lam(maxlat,maxlon,&maxX,&maxY,&proj);
        if (fabs(maxX-minX)>0.0001) {
                float tmp;
                if (minX>maxX) {
                        tmp=minX; minX=maxX; maxX=tmp;
                        tmp=minY; minY=maxY; maxY=tmp;
                }
        } else {
                if (minY>maxY) {
                        float tmp;
                        tmp=minX; minX=maxX; maxX=tmp;
                        tmp=minY; minY=maxY; maxY=tmp;
                }
        }
        fprintf(stderr,"%.2f,%.2f -> %.2f,%.2f\n",
                        minX/1000.,minY/1000.,maxX/1000.,maxY/1000.);
        if (fabs(maxX-minX)>0.0001) {
                double A;
                A=atan2(maxY-minY,maxX-minX);
                DX=hinc*cos(A);
                DY=hinc*sin(A);
        } else { DY=hinc; DX=0.0; }
        fprintf(stderr,"DX=%.4f DY=%.4f\n",DX,DY);
        minX/=1000.;
        maxX/=1000.;
        minY/=1000.;
        maxY/=1000.;

        Nit=sqrt((maxX-minX)*(maxX-minX)+((maxY-minY)*(maxY-minY)))/hinc;

        lay=0;
        fscanf(fgeo,"%lf %lf %f %lf %lf %s %lf %lf\n",
                        &lon,&lat,&Z,&pert,&sd,slay,&dhit,&res);
        while (!feof(fgeo)) {
                int l;
                pt=(struct pt_t*)realloc(pt,(Nblk+1)*sizeof(struct pt_t));
                sscanf(slay+5,"%d",&l);
                if (l>lay) {
                        lay=l;
                }
                geo2lam(lat,lon,&X,&Y,&proj);
                pt[Nblk].pert=pert;
                pt[Nblk].sd=sd;
                pt[Nblk].lay=lay;
                pt[Nblk].hit=dhit;
                pt[Nblk].res=res;
                pt[Nblk].X=X/1000.;
                pt[Nblk].Y=Y/1000.;
                pt[Nblk].Z=Z;
                fscanf(fgeo,"%lf %lf %f %lf %lf %s %lf %lf\n",
                                &lon,&lat,&Z,&pert,&sd,slay,&dhit,&res);
                Nblk++;
        }
        fclose(fgeo);
        fprintf(stderr,"Nblk=%d\nNit=%d\nNlay=%d\n",Nblk,Nit,lay);

        sprintf(geoname,"%s-c.geo",basename);
        fprintf(stderr,"Output: %s\n",geoname);
        if ((fgeo=fopen(geoname,"wt"))==NULL) {
                perror(geoname);
                exit(0);
        }
        lay++;

        depth=0.;
        fscanf(flay,"%*d %lf %*d %lf %lf %*f %*f",&X,&Y,&H);
        depth+=H;
        z=depth;
        fprintf(stderr,"Layer: ");
        for (i=1; i<lay; i++) {
                double dist2,z;
                int startlay=-1,stoplay=-1,k;
                for (j=0; j<Nblk; j++) {
                        if (i==pt[j].lay)
                                if (startlay<0) startlay=j;
                        if (i!=pt[j].lay)
                                if (startlay>=0) { stoplay=j; break; }
                }
                if ((stoplay<0)&&(startlay>=0)) stoplay=Nblk;
                fscanf(flay,"%*d %lf %*d %lf %lf %*f %*f",&X,&Y,&H);
                width=0.5*wfact*(X+Y);
                fprintf(stderr,"%d(%3d,%3d) [%3.1fkm] Z=%.1f..%.1f\n       ",
                                i,startlay,stoplay,width,depth,depth+H);
                for (L=0.,X=minX,Y=minY,k=0; k<Nit; k++,L+=hinc,X+=DX,Y+=DY) {
                        dhit=sum=sd=res=pert=0.;
                        for (j=startlay; j<stoplay; j++) {
                                dist2=(X-pt[j].X)*(X-pt[j].X)+(Y-pt[j].Y)*(Y-pt[j].Y);
                                Z=pt[j].Z;
                                if (dist2>2.*width*width) continue;
                                dist2=exp(-4.0*dist2/(width*width));
                                sum+=dist2;
                                pert+=dist2*pt[j].pert;
                                sd  +=dist2*pt[j].sd  ;
                                dhit +=dist2*pt[j].hit ;
                                res +=dist2*pt[j].res ;
                        }
                        if (dhit<5.0) continue;
                        dhit/=sum;
                        sd/=sum;
                        pert/=sum;
                        res/=sum;
                        for (z=depth; z<depth+H; z+=vinc)
                                fprintf(fgeo,"%10e %10e %10e %10e %10e layer%d %10f %10f\n",
                                                L,0.0,z,pert,sd,i,dhit,res);
                }
                depth+=H;
        }
        fprintf(stderr,"\nEnd point at X=%.2f Y=%.2f L=%.2f\n",X-minX,Y-minY,L);
        fclose(fgeo);
        fclose(flay);
        return(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s -n name -i hinc -z vinc -w width -f lat1 lon1 lat2 lon2\n",argv[0]);
        return(1);
}