dat2matrix.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
 *
 */
/* trace de rai dans un modele "bloc" et eciture du fichier
 * system-name (binaire)
 * */

/* revblk[numero de bloc vrai]=indice de G */
/* numblk[indice de G]=numero de bloc vrai */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <eve.h>
#include <sismoutil.h>
#include "tomo3d.h"

#define RE 6370.

int main(int argc, char **argv)
{
        struct block_t * blk=NULL;
        struct event_t *eve=NULL;
        struct event_t *cur=NULL;
        struct ray_t *ray=NULL;
        struct ray_t *baseray=NULL;
        float **G;
        float *d,*sd;
        FILE *f;
        int MINHIT,MINRAY=3;
        double DS=0.5;
        int Nblk,nbb,num,lastnum,ip,Nray=0,OK;
        char *basename=NULL,*resname=NULL;
        char stacode[1000][5];
        char sysname[255],numname[255],blkname[255],hitname[255];
        int Nsta=0,is,layer1=0,*numblk,*revblk;
        int j,i,stanum;
        int Neve=0,*hit=NULL;
        double dx,dy,dz,X,Y;
        float sigma2=0.0;
        struct lambert proj;
        int norm=0;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit; }
        initproj(PROJ_UTM,-2.5,0.0,0.0,&proj);
        if (argc<2) goto usage;
        for (i=1; i<argc; i++) {
                if (argv[i][0]=='-') {
                        switch (argv[i][1]) {
                                case 'N' :
                                        norm=1;
                                        break;
                                case 'i' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%lf",&DS)) goto error;
                                        break;
                                case 'e' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%d",&MINRAY)) goto error;
                                        break;
                                case 'm' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%d",&MINHIT)) goto error;
                                        break;
                                case 's' :
                                        if (argc<=i+1) goto error;
                                        if (!sscanf(argv[++i],"%f",&sigma2)) goto error;
                                        break;
                                case 'n' :
                                        if (argc<=i+1) goto error;
                                        basename=argv[++i];
                                        break;
                                case 'r' :
                                        if (argc<=i+1) goto error;
                                        resname=argv[++i];
                                        break;
                                default  : goto error;
                        }
                } else goto error;
        }
        if (!basename) goto usage;
        if (!resname) goto usage;

        sprintf(hitname,"%s.hit",basename);
        sprintf(sysname,"%s.sys",basename);
        sprintf(numname,"%s.num",basename);
        sprintf(blkname,"%s.blk",basename);
        fprintf(stderr,"MINHIT=%d\nDS=%.3f\ns2=%.3f\nMINRAY=%d\n",MINHIT,DS,sigma2,MINRAY);
        if ((blk=readblock(blkname, &nbb))==NULL) exit(0);
        hit=(int *)calloc(nbb,sizeof(int));
        for (layer1=0,X=blk[0].z; layer1<nbb; layer1++) {
                if (blk[layer1].z>X+DS) break;
                X=blk[layer1].z;
        }
        fprintf(stderr,"Layer 1 ending at %d\n",layer1);
        fprintf(stderr,"Max depth=%.1f\n",blk[nbb-1].z+blk[nbb-1].dz);
        fprintf(stderr,"Reading %s...",resname);
        fflush(stderr);
        if ((eve=readevent(resname))==NULL) exit(0);

        /* normalisation */
        if (norm) {
                fprintf(stderr," Normalisation...\n");
                for (cur=eve; cur; cur=cur->next) {
                        float locresref=0.;
                        for (i=0; i<cur->Npick; i++) locresref+=cur->pick[i].res;
                        locresref/=(float)cur->Npick;
                        for (i=0; i<cur->Npick; i++)
                                cur->pick[i].res-=locresref;
                }
        }

        fprintf(stderr,"OK\n");

        fprintf(stderr,"Tracing rays ...");
        for (cur=eve; cur; cur=cur->next) {
                if (cur->Npick<MINHIT) continue;
                for (OK=0,ip=0; ip<cur->Npick; ip++) {
                        stanum=-1;
                        for (is=0; is<Nsta; is++) {
                                if (!strcmp(cur->pick[ip].sta->code,stacode[is])) {
                                        stanum=is;
                                        break;
                                }
                        }
                        if (stanum<0) {
                                stanum=Nsta;
                                strcpy(stacode[Nsta],cur->pick[ip].sta->code);
                                Nsta++;
                        }
                        geo2lam(cur->pick[ip].sta->lat,cur->pick[ip].sta->lon,&X,&Y, &proj);
                        lastnum=num=blocknum(X/1000.,Y/1000.,0.0,blk,nbb,0);
                        if (num<0) {       /* is the station in the target volume? */
                                fprintf(stderr,"starting out!!\n");
                                continue;
                        }
                        if (num<layer1) num=stanum;
                        OK=1;
                        /* yes: new ray */
                        hit[num]++;
                        baseray=(struct ray_t *)realloc(baseray,
                                        (Nray+1)*sizeof(struct ray_t));
                        ray=&(baseray[Nray]);
                        Nray++;
                        ray->stanum=stanum;
                        ray->block=num;
                        /* m -> km */
                        ray->x=X/1000.;
                        ray->y=Y/1000.;
                        ray->z=0.001;
                        /*
                         * fprintf(stderr,
                         * "Source %4.4s X=%.0f Y=%.0f: lat=%5.1f lon=%5.1f baz=%5.1f\n",
                         * cur->pick[ip].sta->code,
                         * ray->x,ray->y,
                         * cur->pick[ip].sta->lat,
                         * cur->pick[ip].sta->lon,
                         * cur->pick[ip].baz);
                         * */
                        ray->evenum=Neve;
                        ray->i=(float)asin((blk[num].V
                                                *(cur->pick[ip].slow)*180.)/(M_PI*RE));
                        ray->p=sin(ray->i)/(blk[num].V);
                        ray->az=M_PI*(cur->pick[ip].baz)/180.;
                        ray->L=0.;
                        ray->res=0.0;
                        ray->next=NULL;
                        while ((num=blocknum(ray->x,ray->y,ray->z,blk,nbb,lastnum))>0) {
                                if (num<layer1) num=stanum;
                                ray->i=asin(ray->p*(blk[num].V));
                                dx=DS*sin(ray->az)*sin(ray->i);
                                dy=DS*cos(ray->az)*sin(ray->i);
                                dz=DS*cos(ray->i);
                                if (dx>3.0*DS) fprintf(stderr,"!!! dx=%f\n",dx);
                                if (dy>3.0*DS) fprintf(stderr,"!!! dy=%f\n",dy);
                                ray->x+=dx;
                                ray->y+=dy;
                                ray->z+=dz;
                                ray->block=num;
                                baseray[Nray-1].res+=ray->L;
                                ray->L+=DS/(blk[num].V);
                                if (num!=lastnum) {
                                        /*
                                         * fprintf(stderr,"Changement de bloc...%d (%3.1f)->%d\n",
                                         * lastnum,ray->L,num);
                                         * */
                                        hit[num]++;
                                        ray->next=(struct ray_t *)malloc(sizeof(struct ray_t));
                                        memcpy(ray->next,ray,sizeof(struct ray_t));
                                        ray->next->next=NULL;
                                        ray=ray->next;
                                        ray->L=0.0;
                                        lastnum=num;
                                }
                        }
                        baseray[Nray-1].res=cur->pick[ip].res;
                        if (ray->z<blk[nbb-1].z+blk[nbb-1].dz) {
                                fprintf(stderr,"Zout<Zmax Nray=%d\n",Nray-1);
                                fprintf(stderr,"x=%f y=%f z=%f\n",ray->x,ray->y,ray->z);
                                exit(0);
                        }
                        /*
                         * fprintf(stderr,"Out point: %f %f %f\n",ray->x,ray->y,ray->z);
                         * */
                }
                if (OK) Neve++;
        }

        fprintf(stderr,"OK\n");
        fprintf(stderr,"%d rays, %d stations, %d events\n",Nray, Nsta, Neve);
        fprintf(stderr,"Allocation..."); fflush(stderr);
        numblk=(int*)calloc(nbb,sizeof(int));
        revblk=(int*)calloc(nbb,sizeof(int));
        d=(float*)calloc(Nray,sizeof(float));
        sd=(float*)calloc(Nray,sizeof(float));
        for (i=0; i<Nray; i++) sd[i]=sigma2;
        fprintf(stderr,"OK\n");
        fprintf(stderr,"OK\nLooking for visited blocks..."); fflush(stderr);
        for (Nblk=0,j=0; j<nbb; j++)
                for (i=0; i<Nray; i++)
                        if (hit[j]>=MINHIT) {
                                revblk[j]=Nblk;
                                numblk[Nblk++]=j;
                                break;
                        }
        fprintf(stderr,"%d\n",Nblk);

        fprintf(stderr,"Init..."); fflush(stderr);
        G=(float**)calloc(Nray,sizeof(float*));
        for (i=0; i<Nray; i++) {
                if ((G[i]=(float*)calloc(Nblk,sizeof(float)))==NULL) {
                        fprintf(stderr,"malloc error\n");
                        exit(1);
                }
        }
        for (i=0; i<Nray; i++) {
                /*d[i]=baseray[i].res;*/
                /* on range dans d[i] le temps total passe dans le modele */
                for (ray=&(baseray[i]); ray; ray=ray->next) d[i]+=ray->L;

                for (ray=&(baseray[i]); ray; ray=ray->next)
                        G[i][revblk[ray->block]]=ray->L/d[i];
                sd[i]/=d[i];

                /* d[i] est en fait le residu raporte au temps total */
                d[i]=baseray[i].res/d[i];
        }

        fprintf(stderr,"Writing..."); fflush(stderr);
        if (NULL==(f=fopen(sysname,"wb"))) {
                perror(sysname);
                exit(1);
        }
        fprintf(stderr,"[header]... "); fflush(stderr);
        fwrite(&Nray,sizeof(int),1,f);
        fwrite(&Nblk,sizeof(int),1,f);
        fprintf(stderr,"[G]... "); fflush(stderr);
        for (i=0; i<Nray; i++)
                for (j=0; j<Nblk; j++)
                        fwrite(&(G[i][j]),sizeof(float),1,f);
        fprintf(stderr,"[d]... "); fflush(stderr);
        fwrite(d,sizeof(float),Nray,f);
        fwrite(sd,sizeof(float),Nray,f);
        fprintf(stderr,"[damp]... "); fflush(stderr);
        layer1=j=0;
        X=blk[0].z;
        for (i=0; i<Nblk; i++)
                fwrite(&(blk[numblk[i]].damp),sizeof(float),1,f);
        fclose(f);
        fprintf(stderr,"[numblk]... "); fflush(stderr);
        if (NULL==(f=fopen(numname,"wb"))) {
                perror(numname);
                exit(1);
        }
        for (i=0; i<Nblk; i++) fprintf(f,"%d\n",numblk[i]);
        fclose(f);
        fprintf(stderr,"[hitblk]... "); fflush(stderr);
        if (NULL==(f=fopen(hitname,"wt"))) {
                perror(hitname);
                exit(1);
        }
        for (i=0; i<Nblk; i++) fprintf(f,"%d\n",hit[numblk[i]]);
        fclose(f);
        fprintf(stderr,"Done\n");
        return(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s -n name -r resfile [-N] [-i inc (0.5)] [-e minray] [-m minhit (3)]\n",
                        argv[0]);
        return(1);
}