ray2nmatrix.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
 *
 */
/* revblk[numero de noeud vrai]=indice de G */
/* numblk[indice de G]=numero de noeud vrai */

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

#define RE 6370.
#define OK fprintf(stderr,"FILE %s, LINE %d\n",__FILE__,__LINE__);
#define EPS 1.e-8

int COUNT=0,LOOPS=0;

struct bl_t {
        int num;
        float T;
        float t2;
        struct bl_t *next;
};

int main(int argc, char **argv)
{
        struct block_t * blk=NULL;
        struct path_t * path=NULL;
        struct bl_t *blt=NULL;
        float **G;
        float *d,*sd;
        double FZ=8.0*8.0;
        float *X=NULL,*Y=NULL,*Z=NULL,*T=NULL;
        FILE *f,*fnum;
        float DS,maxdepth;
        int MINHIT=3,Nray,i,j,nbb,layer1,NAFF;
        int NSEGMAX;
        int *hit=NULL,*numblk=NULL,*revblk=NULL;
        int num,Nblk=0;
        float sigma2=0.1;
        char sysname[255],numname[255],blkname[255];
        char hitname[255],rayname[255];
        char *basename=NULL;
        struct path_t *p;
        struct bl_t *bltmp;
        int *useblk=NULL;
        int Nused=0;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit(0); }
        if (argc<2) goto usage;
        for (i=1; i<argc; i++) {
                if (argv[i][0]=='-') {
                        switch (argv[i][1]) {
                                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;
                                default  :
                                        goto error;
                        }
                } else goto error;
        }
        if (!basename) goto usage;

        sprintf(hitname,"%s.hit",basename);
        sprintf(sysname,"%s.sys",basename);
        sprintf(numname,"%s.num",basename);
        sprintf(blkname,"%s.blk",basename);
        sprintf(rayname,"%s.ray",basename);

        if ((blk=readblock(blkname, &nbb))==NULL) exit(1);
        for (layer1=0,DS=blk[0].z; layer1<nbb; layer1++) {
                if (blk[layer1].z>DS+blk[0].z/10.) break;
                DS=blk[layer1].z;
        }
        useblk=(int*)calloc(nbb,sizeof(int));
        fprintf(stderr,"Reading %s...",numname);
        if ((fnum=fopen(numname,"rt"))==NULL) {
                perror(numname);
                exit(0);
        }
        fscanf(fnum,"%d",&(useblk[Nused++]));
        while (!feof(fnum)) if (!fscanf(fnum,"%d",&(useblk[Nused++]))) break;
        fclose(fnum);
        Nused--;
        fprintf(stderr,"used nodes: %d\n",Nused);
        fprintf(stderr,"Layer 1 ending at %d, ",layer1);
        fprintf(stderr,"max depth=%.1f\n",maxdepth=blk[nbb-1].z+blk[nbb-1].dz);

        fprintf(stderr,"Reading %s...",rayname);
        if ((f=fopen(rayname,"rb"))==NULL) {
                perror(rayname);
                exit(1);
        }
        fread(&Nray, sizeof(int), 1, f);
        fread(&DS, sizeof(float), 1, f);
        /*path=readpath(rayname,&Nray,&DS);*/
        NAFF=Nray/10;
        if (!NAFF) NAFF=10;
        fprintf(stderr,"%d rays, DS=%f\n",Nray,DS);

        fprintf(stderr,"Allocation..."); fflush(stderr);
        hit   =(int*)calloc(nbb,sizeof(int));
        numblk=(int*)calloc(nbb,sizeof(int));
        revblk=(int*)calloc(nbb,sizeof(int));
        d= (float*)calloc(Nray,sizeof(float));
        sd=(float*)calloc(Nray,sizeof(float));
        blt=(struct bl_t*)malloc(Nray*sizeof(struct bl_t));
        if ((!hit)||(!numblk)||(!revblk)||(!d)||(!sd)||(!blt)) {
                fprintf(stderr,"alloc error\n");
                exit(1);
        }
        for (i=0; i<nbb; i++) revblk[i]=numblk[i]=-1;

        fprintf(stderr,"OK\nBlock pass..."); fflush(stderr);
        /* allocation for all the ray headers (= path) */
        path=(struct path_t *)malloc(Nray*sizeof(struct path_t));
        NSEGMAX=-1;
        ptime(0,1);
        for (i=0; i<Nray; i++) {
                double dist=0.0;
                double Tdist=0.0;
                int k,kk;
                p=&(path[i]);
                /* read the ray header */
                fread(p,sizeof(struct path_t),1,f);
                if (p->nseg>NSEGMAX) {
                        fprintf(stderr," +");
                        NSEGMAX=p->nseg+50;
                        X=realloc(X,sizeof(float)*NSEGMAX);
                        Y=realloc(Y,sizeof(float)*NSEGMAX);
                        Z=realloc(Z,sizeof(float)*NSEGMAX);
                        T=realloc(T,sizeof(float)*NSEGMAX);
                        if ((!X)||(!Y)||(!Z)||(!T)) {
                                fprintf(stderr,"alloc error\n");
                                exit(1);
                        }
                }

                /* read the ray path and time along segments DS */
                fread(X,sizeof(float),p->nseg,f);
                fread(Y,sizeof(float),p->nseg,f);
                fread(Z,sizeof(float),p->nseg,f);
                fread(T,sizeof(float),p->nseg,f);

                if (!(i%NAFF)) fprintf(stderr," %d",i);

                /* blt[i] = starting block for ray i = block stanum */
                bltmp=&(blt[i]);
                bltmp->num=num=p->stanum;
                bltmp->T=0.0;
                bltmp->t2=0.0;
                bltmp->next=NULL;
                hit[num]++;
                /* the ray is going back from the station towards the source */
                for (j=0; j<p->nseg; j++,LOOPS++) {
                        /* below last layer ? */
                        if (Z[j]>maxdepth) break;
                        /* in the first layer ? */
                        if (Z[j]<=blk[0].dz) {
                                bltmp->T+=T[j];
                                continue;
                        }
                        COUNT++;
                        if (num<0) {      /* out of the model ? */
                                fprintf(stderr,"\nout: ray %d: %f %f %f\n",
                                                i,X[j],Y[j],Z[j]);
                                exit(1);
                        }

                        Tdist=0.;
                        for (kk=0; kk<Nused; kk++) {
                                float tmp;
                                k=useblk[kk];
                                dist=0.;

                                tmp=(blk[k].z-Z[j]);
                                dist+=tmp*tmp;
                                if (dist>5.0*FZ) continue;
                                tmp=(blk[k].x-X[j]);
                                dist+=tmp*tmp;
                                if (dist>5.0*FZ) continue;
                                tmp=(blk[k].y-Y[j]);
                                dist+=tmp*tmp;
                                if (dist>5.0*FZ) continue;

                                dist=exp(-0.5*dist/FZ);
                                if (dist<EPS) continue;
                                Tdist+=dist;

                                for (bltmp=&(blt[i]); bltmp->next; bltmp=bltmp->next)
                                        if (bltmp->num==k) break;
                                if (bltmp->num!=k) {
                                        bltmp->next=(struct bl_t*)calloc(1,sizeof(struct bl_t));
                                        if (!bltmp->next) {
                                                fprintf(stderr,"alloc error\n");
                                                exit(1);
                                        }
                                        bltmp=bltmp->next;
                                        bltmp->T=0.;
                                        bltmp->t2=0.;
                                        bltmp->next=NULL;
                                        bltmp->num=k;
                                }
                                bltmp->t2=dist*T[j];
                        }
                        if (Tdist<0.5*EPS) continue;
                        for (k=0,bltmp=blt[i].next; bltmp; bltmp=bltmp->next,k++) {
                                hit[bltmp->num]++;
                                bltmp->T+=bltmp->t2/Tdist;
                                bltmp->t2=0.;
                        }
                }
        }
        ptime(0,2);
        fprintf(stderr,"\nLoops=%d\n",LOOPS);
        fprintf(stderr,"Reordering blocks..."); fflush(stderr);
        for (Nblk=0,j=0; j<nbb; j++)
                if ((hit[j]>=MINHIT)||((hit[j]>0)&&(j<layer1))) {
                        revblk[j]=Nblk;
                        numblk[Nblk]=j;
                        Nblk++;
                }
        fprintf(stderr,"%d\n",Nblk);
        /* The actual number of visited blocks now Nblk */

        fprintf(stderr,"Init..."); fflush(stderr);
        G=(float**)calloc(Nray,sizeof(float*));
        if (!G) {
                fprintf(stderr,"alloc error\n");
                exit(1);
        }
        for (i=0; i<Nray; i++) {
                if (!(i%NAFF)) fprintf(stderr,"%d ",i);
                if ((G[i]=(float*)calloc(Nblk,sizeof(float)))==NULL) {
                        fprintf(stderr,"\nmalloc error\n");
                        exit(1);
                }
        }

        fprintf(stderr,"\nComputing G[i][j]...");
        for (i=0; i<Nray; i++) {
                struct bl_t *bltmp=NULL;

                /* d[i]=temps total passe par le rai i dans le modele */
                d[i]=0.0;
                for (bltmp=&(blt[i]); bltmp; bltmp=bltmp->next) d[i]+=bltmp->T;
                d[i]=1.0;   /* NON ca veut dire que TOUTES les donnees ont
                               le meme poids relatif = division de chaque ligne du
                               systeme par la somme des coeff de la ligne de G */

                /* distribution des perturbations dans les blocs */
                for (bltmp=&(blt[i]); bltmp; bltmp=bltmp->next) {
                        if (revblk[bltmp->num]>=0)
                                G[i][revblk[bltmp->num]]=bltmp->T/d[i];
                }

                /* pareil sur les perturbations et les erreurs */
                sd[i]=sigma2/d[i];
                d[i]=path[i].res/d[i];
        }

        ptime(0,1);
        fprintf(stderr,"\nWriting...[numblk]... "); fflush(stderr);
        if (NULL==(f=fopen(numname,"wt"))) {
                perror(numname);
                exit(1);
        }
        for (i=0; i<Nblk; i++) fprintf(f,"%d\n",numblk[i]);
        fclose(f);
        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);
                fflush(f);
        }
        fprintf(stderr,"[d]... "); fflush(stderr);
        fwrite(d,sizeof(float),Nray,f);
        fprintf(stderr,"[sd]... "); fflush(stderr);
        fwrite(sd,sizeof(float),Nray,f);
        fprintf(stderr,"[damp]... "); fflush(stderr);
        for (i=0; i<Nblk; i++)
                fwrite(&(blk[numblk[i]].damp),sizeof(float),1,f);
        fclose(f);
        fprintf(stderr,"\nWriting..."); fflush(stderr);
        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);
        ptime(0,2);
        fprintf(stderr,"\n");
        return(0);
error:
        fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
        fprintf(stderr,"%s -n name [-s sigma2] [-m minhit (3)]\n",
                        argv[0]);
        return(1);
}