postres.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 <eve.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <sismoutil.h>
#include "tomo3d.h"

double homo(double x)
{ return(8.0); }

int main(int argc, char **argv)
{
        FILE *f;
        char rayname[255],outname[255];
        char *basename=NULL;
        struct event_t *cur=NULL;
        struct data_t *data=NULL;
        struct data_t *dcur=NULL;
        float *res=NULL;
        int Nray=0,i,MINRAY=3,ires=0;
        char *refmodname=NULL;
        float P2Sconv;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit(0); }
        init_mod(NULL,0.);
        if (argc<2) goto usage;
        for (i=1; i<argc; i++) {
                if (argv[i][0]=='-') {
                        switch (argv[i][1]) {
                                case 'c' :
                                        if (argc<=i+2) goto error;
                                        refmodname=argv[++i];
                                        if (!sscanf(argv[++i],"%f",&P2Sconv)) goto error;
                                        init_mod(refmodname,P2Sconv);
                                        break;
                                case 'm' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%d",&MINRAY)) goto error;
                                        break;
                                case 's' :
                                case 'h' :
                                case 'p' :
                                        if (argc<=i+1) goto error;
                                        if (data) {
                                                data->next=(struct data_t*)malloc(sizeof(struct data_t));
                                                dcur=data->next;
                                        } else dcur=data=(struct data_t*)malloc(sizeof(struct data_t));
                                        dcur->resname=argv[i+1];
                                        switch (argv[i][1]) {
                                                case 'h' :
                                                        dcur->vel=homo;
                                                        fprintf(stderr,"Pdata: %s\n",dcur->resname);
                                                        break;
                                                case 'p' :
                                                        dcur->vel=Pvelocity;
                                                        fprintf(stderr,"Pdata (IASP): %s\n",dcur->resname);
                                                        break;
                                                case 's' :
                                                        dcur->vel=Svelocity;
                                                        fprintf(stderr,"Sdata (IASP): %s\n",dcur->resname);
                                                        break;
                                        }
                                        dcur->next=NULL;
                                        dcur->eve=NULL;
                                        i++;
                                        break;
                                case 'n' :
                                        if (argc<=i+1) goto error;
                                        basename=argv[++i];
                                        break;
                                default  :
                                        goto error;
                        }
                } else goto error;
        }
        if (!basename) goto usage;

        fprintf(stderr,"MINRAY=%d\n",MINRAY);

        sprintf(outname,"%s.out",basename);
        sprintf(rayname,"%s.ray",basename);


        for (dcur=data; dcur; dcur=dcur->next) {
                fprintf(stderr,"Reading %s...",dcur->resname);
                if ((dcur->eve=readevent(dcur->resname))==NULL) exit(0);
        }
        fprintf(stderr,"Selecting events...\n");
        Nray=0;
        for (dcur=data; dcur; dcur=dcur->next)
                for (cur=dcur->eve; cur; cur=cur->next) {
                        if (cur->Npick<MINRAY) cur->Npick=0; else Nray+=cur->Npick;
                }
        fprintf(stderr,"Cleaning/numbering events...\n");
        for (dcur=data; dcur; dcur=dcur->next)
                dcur->eve=cleanevent(dcur->eve);

        fprintf(stderr,"Estimated Nray=%d\n",Nray);

        fprintf(stderr,"Reading %s...",rayname);
        if ((f=fopen(rayname,"rb"))==NULL) {
                perror(rayname);
                exit(1);
        }
        fread(&i, sizeof(int), 1, f);
        fclose(f);
        fprintf(stderr,"\nEffective Nray=%d\n",i);
        if (i!=Nray) exit(0);

        fprintf(stderr,"Reading %s...",outname);
        if ((f=fopen(outname,"rb"))==NULL) {
                perror(outname);
                exit(1);
        }
        res= (float*)calloc(Nray,sizeof(float));
        if (!res) {
                fprintf(stderr,"alloc error\n");
                exit(1);
        }
        for (i=0; i<Nray; i++)
                if (!fscanf(f,"%*d %*f %f",&(res[i]))) {
                        fprintf(stderr,"error reading %s (i=%d)\n",outname,i);
                        exit(0);
                }
        fclose(f);
        fprintf(stderr,"\n");

        ires=0;
        for (dcur=data; dcur; dcur=dcur->next) {
                for (cur=dcur->eve; cur; cur=cur->next) {
                        if (cur->Npick<MINRAY) continue;
                        for (i=0; i<cur->Npick; ires++,i++) {
                                cur->pick[i].res=res[ires];
                                cur->pick[i].obs=cur->pick[i].cal+cur->pick[i].res;
                                strcpy(cur->eve->name,"residual");
                        }
                }
        }

        for (dcur=data; dcur; dcur=dcur->next) {
                fprintf(stdout,"RESIDUAL %s\n",dcur->resname);
                writedat(stdout,dcur->eve);
        }

        return(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s -n name [-c refmod-file P2Sconv] [-p Presfile] [-p ...] [-s Sresfile] [-s ...] [-m minray]\n",
                        argv[0]);
        return(1);
}