ray2bmatrixnew.c

射线追踪程序

/* 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.

extern int SAVED;
int COUNT=0,LOOPS=0;

struct bl_t {
        int num,evenum;
        float T;
        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;
        float *X=NULL,*Y=NULL,*Z=NULL,*T=NULL;
        float **sG=NULL;
        int *nG=NULL;
        int ACH=0;
        int NORM=0;
        FILE *f;
        float DS,maxdepth;
        int MINHIT=3,Nray,i,j,nbb,layer1,NAFF;
        int NSEGMAX;
        int *hit=NULL,*numblk=NULL,*revblk=NULL;
        int num,lastblk,Nblk=0;
        float sigma2=0.1;
        char sysname[255],numname[255],blkname[255];
        char hitname[255],rayname[255];
        char *basename=NULL;
        int cache=0;
        struct path_t *p;
        struct bl_t *bltmp;
        int Neve=0;
        argv=parseoptions(&argc,argv);
        if (argc<0) { argc*=-1; exit; }
        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 'A' :
                                        ACH=1;
                                        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;
                                default  :
                                        goto error;
                        }
                } else goto error;
        }
        if ((ACH)&&(NORM)) {
                fprintf(stderr,"ACH xor NORM ???\n");
                exit(1);
        }
        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);
        qsort((void*)blk, nbb, sizeof(struct block_t),blkcmp);
        for (layer1=0,DS=blk[0].z; layer1<nbb; layer1++) {
                if (blk[layer1].z>DS+blk[0].z/10.) break;
                DS=blk[layer1].z;
        }
        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));
        SAVED=0;
        NSEGMAX=-1;
        ptime(0,1);
        for (i=0; i<Nray; i++) {
                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=lastblk=p->stanum;
                bltmp->evenum=p->evenum;
                bltmp->T=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++) {
                        /* in which block is the ray ? */

                        if (Z[j]>maxdepth) break; /* below last layer ? */

                        if (Z[j]<=blk[0].dz) num=p->stanum; /* in the first layer ? */
                        else cache=num=blocknum(X[j],Y[j],Z[j],blk,nbb,cache);
                        COUNT++;
                        if (num<0) {      /* out of the model ? */
                                fprintf(stderr,"\n\n\n\n\n\n\n\nOUT ray %d: %f %f %f\n\n\n\n\n",
                                                i,X[j],Y[j],Z[j]);
                                exit(1);
                        }
                        if (lastblk!=num) {
                                /* the ray goes into an other block : one more */
                                if (!bltmp->next) {
                                        int inbl=0;
                                        struct bl_t *tmpbl=bltmp;
                                        for (inbl=0; inbl<20; inbl++) {
                                                tmpbl->next=(struct bl_t*)calloc(1,sizeof(struct bl_t));
                                                tmpbl=tmpbl->next;
                                        }
                                        if (!bltmp->next) {
                                                fprintf(stderr,"alloc error\n");
                                                exit(1);
                                        }
                                }
                                bltmp=bltmp->next;
                                bltmp->next=NULL;
                                bltmp->num=lastblk=num;
                                bltmp->evenum=p->evenum;
                                hit[num]++;
                                bltmp->T=0.;
                                /* bltmp is now the actual block associated to the segment */
                        }
                        bltmp->T+=T[j];
                }
        }
        ptime(0,2);
        fprintf(stderr,"\nLoops=%d\n",LOOPS);
        fprintf(stderr,"Saved scans=%d/%d (%.2f%%)\n",
                        COUNT-SAVED,COUNT,100.*(float)(COUNT-SAVED)/(float)COUNT);
        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);
                /* revblk[numero de bloc vrai]=indice de G 
                 * numblk[indice de G]=numero de bloc vrai 
                 * nbb=nb de blocs dans la grille */
                /* The actual number of visited blocks now Nblk */

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

                if (ACH) {
                        int leve=0;
                        Neve=0;
                        blt[0].evenum=0;
                        for (i=0; i<Nray; i++) {
                                if (blt[i].evenum!=leve) {
                                        leve=blt[i].evenum;
                                        Neve++;
                                }
                                blt[i].evenum=Neve;
                        }
                        Neve++;
                        fprintf(stderr,"  ACH: Neve=%d\n",Neve);

                        sG=(float**)calloc(Neve,sizeof(float*));
                        if (!sG) {
                                fprintf(stderr,"sG malloc error\n");
                                exit(1);
                        }
                        for (i=0; i<Neve; i++)
                                if ((sG[i]=(float*)calloc(Nblk,sizeof(float)))==NULL) {
                                        fprintf(stderr,"sG[%d] malloc error\n",i);
                                        exit(1);
                                }

                        nG=(int*)calloc(Neve,sizeof(int));
                        if (!nG) {
                                fprintf(stderr,"nG malloc error\n");
                                exit(1);
                        }
                        for (i=0; i<Nray; i++)  {
                                int leve=blt[i].evenum;
                                if (i) if (leve<blt[i-1].evenum)
                                        fprintf(stderr,"oups i=%d leve=%d\n",i,leve);
                                for (bltmp=&(blt[i]); bltmp; bltmp=bltmp->next) {
                                        int nb;
                                        nb=revblk[bltmp->num];
                                        if (nb<0) continue;
                                        sG[leve][nb]+=bltmp->T;
                                }
                                nG[leve]++;
                        }
                }

                for (i=0; i<Nray; i++) {
                        int leve=blt[i].evenum;
                        struct bl_t *bltmp=NULL;

                        /* d[i]=temps total passe par le rai i dans le modele */
                        if (NORM) {
                                d[i]=0.0;
                                for (bltmp=&(blt[i]); bltmp; bltmp=bltmp->next) d[i]+=bltmp->T;
                        }
                        else d[i]=1.0;

                        /* distribution des perturbations dans les blocs */
                        for (bltmp=&(blt[i]); bltmp; bltmp=bltmp->next) {
                                int nb;
                                nb=revblk[bltmp->num];
                                if (nb<0) continue;
                                G[i][nb]=bltmp->T/d[i];
                                if (ACH) G[i][nb]-=(sG[leve][nb]/(float)nG[leve]);
                        }

                        /* 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..."); 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);
                        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...[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);
                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 [-A] [-s sigma2] [-m minhit (3)]\n",
                                argv[0]);
                return(1);
        }
}