tomoresample.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 <stdio.h>
#include <stdlib.h>

#define UNIT 10000

struct data_t {
        float lat,lon,dep;
        double X,Y;
        int layer;
        float dv,res,hit,err;
};

int main(int argc,char **argv)
{
        struct lambert proj;
        float MINDIST;
        struct data_t *data;
        int Ndata=0,i,ntemp;
        float latmin,latmax,lonmin,lonmax,depmin,depmax;
        float dx=5.0,dx111;
        float lat,lon,dep,dv,res,hit,err;
        char lay[200];
        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 'x' :
                                        sscanf(argv[i+1],"%f",&dx);
                                        break;
                        }
                }
        }
        dx111=dx/111.19;
        MINDIST=4.0*dx;
        initproj(13,46.0,4.0,0.0,&proj);
        ntemp=UNIT;
        data=(struct data_t*)malloc(ntemp*sizeof(struct data_t));
        if (!data) {
                mallocerror();
                exit(1);
        }
        while (8==fscanf(stdin,"%f %f %f %f %f %s %f %f",
                                &lat,&lon,&dep,&dv,&res,lay,&hit,&err)) {
                data[Ndata].lat=lat;
                data[Ndata].lon=lon;
                data[Ndata].dep=dep;
                data[Ndata].dv=dv;
                data[Ndata].res=res;
                data[Ndata].hit=hit;
                data[Ndata].err=err;
                sscanf(lay,"layer%d",&(data[Ndata].layer));
                geo2lam(lat,lon,&(data[Ndata].X),&(data[Ndata].Y),&proj);
                data[Ndata].X/=1000.0;
                data[Ndata].Y/=1000.0;
                Ndata++;
                if (Ndata>=ntemp) {
                        ntemp+=UNIT;
                        data=(struct data_t*)realloc(data,ntemp*sizeof(struct data_t));
                        printlog("%d read\n",Ndata);
                        if (!data) {
                                mallocerror();
                                exit(1);
                        }
                }
        }
        printlog("Ndata=%d\n",Ndata);
        latmax=latmin=data[0].lat;
        lonmax=lonmin=data[0].lon;
        depmax=depmin=data[0].dep;
        for (i=0; i<Ndata; i++) {
                if (data[i].lat<latmin) {
                        latmin=data[i].lat;
                }
                if (data[i].lat>latmax) {
                        latmax=data[i].lat;
                }
                if (data[i].lon<lonmin) {
                        lonmin=data[i].lon;
                }
                if (data[i].lon>lonmax) {
                        lonmax=data[i].lon;
                }
                if (data[i].dep<depmin) {
                        depmin=data[i].dep;
                }
                if (data[i].dep>depmax) {
                        depmax=data[i].dep;
                }
        }
        printlog("lat: %.2f -> %.2f\n",latmin,latmax);
        printlog("lon: %.2f -> %.2f\n",lonmin,lonmax);
        printlog("dep: %.2f -> %.2f\n",depmin,depmax);
        latmin=floor(latmin);
        lonmin=floor(lonmin);
        depmin=0.0;
        latmax=ceil(latmax);
        lonmax=ceil(lonmax);
        depmax=ceil(depmax);
        printlog("LAT: %.2f -> %.2f\n",latmin,latmax);
        printlog("LON: %.2f -> %.2f\n",lonmin,lonmax);
        printlog("DEP: %.2f -> %.2f\n",depmin,depmax);
        for (lat=latmin; lat<latmax; lat+=dx111) {
                for (lon=lonmin; lon<lonmax; lon+=dx111) {
                        for (dep=depmin; dep<depmax; dep+=dx) {
                                double P,X,Y,sum,d1,dist;
                                geo2lam(lat,lon,&X,&Y,&proj);
                                X/=1000.0;
                                Y/=1000.0;
                                dv=err=hit=res=sum=0.0;
                                for (i=0; i<Ndata; i++) {
                                        dist=0.0;

                                        d1=(X-data[i].X);
                                        dist+=d1*d1;
                                        if (dist>MINDIST) {
                                                continue;
                                        }

                                        d1=(Y-data[i].Y);
                                        dist+=d1*d1;
                                        if (dist>MINDIST) {
                                                continue;
                                        }

                                        d1+=(dep-data[i].dep);
                                        dist+=d1*d1;
                                        if (dist>MINDIST) {
                                                continue;
                                        }

                                        dist/=dx;
                                        dist/=dx;
                                        P=exp(-4.0*dist);
                                        dv+=P*data[i].dv;
                                        err+=P*data[i].err;
                                        res+=P*data[i].res;
                                        hit+=P*data[i].hit;
                                        sum+=P;
                                }
                                if (sum>1.0e-8) {
                                        printf("%f %f %f %f\n",dv/sum,res/sum,hit/sum,err/sum);
                                } else {
                                        printf("-100.0 -100.0 -100.0 -100.0\n");
                                }
                        }
                        printlog("done: %f %f %f\n",lat,lon,dep);
                }
        }

}