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

/*#undef NTHREADS*/

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

clock_t __lclock=0;
#define PROGRESS(k,N) { if (clock()>__lclock+CLOCKS_PER_SEC) { fprintf(stderr,"%6ld/%6ld",k,N); __lclock=clock(); } }

#define RE 6370.
#define OK fprintf(stderr,"FILE %s, LINE %d\n",__FILE__,__LINE__);
#define MAXBLR 15

struct bl_t {
        int num,evenum;
        float T;
        struct bl_t *next;
};

#ifdef NTHREADS
#include <pthread.h>

struct pt_param {
        int id;
        struct bl_t *bltmp;
        struct block_t * blk;
        int nbb,i;
        float *X,*Y,*Z,*T;
        float maxdepth;
        struct path_t *p;
};

struct pt_param *pt;
pthread_t *threads = NULL;
int nbt,nt;
int *status;
#endif

void processray(bltmp,blk,nbb,X,Y,Z,T,p,i,maxdepth)
struct bl_t *bltmp;
struct block_t * blk;
int nbb,i;
float *X,*Y,*Z,*T;
float maxdepth;
struct path_t *p;
{
        int j,cache=0,num,lastblk;
        int Nblr=0;
        struct bl_t *reserve=NULL;
        lastblk=bltmp->num;
        for (j=0; j<p->nseg; j++) {

                /* 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);
                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",
                                        i,X[j],Y[j],Z[j]);
                        exit(1);
                }
                if (lastblk!=num) {
                        /* the ray goes into an other block : one more */
                        if (!Nblr) {
                                reserve=(struct bl_t*)malloc(MAXBLR*sizeof(struct bl_t));
                                assert(reserve);
                                Nblr=MAXBLR;
                        }
                        bltmp->next=&(reserve[--Nblr]);
                        bltmp=bltmp->next;
                        bltmp->next=NULL;
                        bltmp->num=lastblk=num;
                        bltmp->evenum=p->evenum;
                        bltmp->T=0.;
                        /* bltmp is now the actual block associated to the segment */
                }
                bltmp->T+=T[j];
        }
}

#ifdef NTHREADS
void * thr_processray(void *pointer)
{
        struct pt_param *p=(struct pt_param*)pointer;
        processray(p->bltmp,p->blk,p->nbb,p->X,p->Y,p->Z,p->T,
                        p->p,p->i,p->maxdepth);
        return(&(p->id));
}
#endif

int main(int argc, char **argv)
{
        struct block_t * blk=NULL;
        struct path_t * path=NULL;
        struct bl_t *blt=NULL;
        struct matrix_t *G;
        float *d,*sd;
#ifndef NTHREADS
        float *X=NULL,*Y=NULL,*Z=NULL,*T=NULL;
#endif
        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],
        hitname[255],rayname[255];
        char *basename=NULL;
        struct path_t *p;
        struct bl_t *bltmp;
        int Neve=0;
        time_t iclock;
        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 '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;

        fprintf(stderr,"=============== ray2Bmatrix ===============\n");

        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=100;
        if (!NAFF) NAFF=10;
        fprintf(stderr,"%d rays @ DS=%.1f\n",Nray,DS);

        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,"Block pass...");
        /* allocation for all the ray headers (= path) */
        path=(struct path_t *)malloc(Nray*sizeof(struct path_t));
        NSEGMAX=-1;