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

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

double compT(int Nseg, float *X, float *Y, float *Z, float *T,
                float ds, float *L,
                float *Ah, float *Av, int nh,
                struct mod3d_t *mod, int modif)
{
        int j,k;
        double Tc=0.,xx,yy,zz,lx,ly,lz;
        double l,H,tmp,t,tl,d;
        double khx,khy,khz,kvz,kvx,kvy;
        H=sqrt((X[Nseg-1]-X[0])*(X[Nseg-1]-X[0])+(Y[Nseg-1]-Y[0])*(Y[Nseg-1]-Y[0]));
        khx=-(Y[Nseg-1]-Y[0])/H;
        khy=+(X[Nseg-1]-X[0])/H;
        khz=0.;
        kvx=-(Z[Nseg-1]/(*L))*((X[Nseg-1]-X[0])/H);
        kvy=-(Z[Nseg-1]/(*L))*((Y[Nseg-1]-Y[0])/H);
        kvz=+H/(*L);
        tl=0.;
        lx=X[0];
        ly=Y[0];
        lz=Z[0];
        for (j=1; j<Nseg; j++) {
                xx=X[j];
                yy=Y[j];
                zz=Z[j];
                l=(double)j*ds*M_PI/(*L);
                for (k=0; k<nh; k++) {
                        tmp=sin(((double)k+0.5)*l)/(double)(k+1);
                        xx+=(khx*Ah[k]+kvx*Av[k])*tmp;
                        yy+=(khy*Ah[k]+kvy*Av[k])*tmp;
                        zz+=(khz*Ah[k]+kvz*Av[k])*tmp;
                }
                d=sqrt((xx-lx)*(xx-lx)+
                                (yy-ly)*(yy-ly)+
                                (zz-lz)*(zz-lz));
                t=d/getvelocity(xx,yy,zz,mod,0);
                tl+=d;
                if (modif) {
                        /*
                         * printf("%f %f %f P\n",xx,yy,zz);
                         * printf("%f %f %f O\n",X[j],Y[j],Z[j]);
                         * */
                        fflush(stdout);
                        X[j]=xx;
                        Y[j]=yy;
                        Z[j]=zz;
                        T[j]=t;
                }
                Tc+=t;
                lx=xx;
                ly=yy;
                lz=zz;
        }
        *L=(float)tl;
        return(Tc);
}

void sorttab(int nparam, float *key, float **tab)
{
        float *p,tmp;
        int i,j;
        for (i=0; i<nparam; i++)
                for (j=i+1; j<nparam; j++)
                        if (key[j]<key[i]) {
                                p=tab[i];
                                tab[i]=tab[j];
                                tab[j]=p;
                                tmp=key[i];
                                key[i]=key[j];
                                key[j]=tmp;
                        }
}

void swtch(float **p1,float **p2,float *v1,float *v2)
{
        float *p,v;
        p=*p1; *p1=*p2; *p2=p;
        v=*v1; *v1=*v2; *v2=v;
}

void findbestSIMPLEX(struct path_t * p,
                float *X, float *Y, float *Z, float *T,int NH,
                struct mod3d_t *mod, float ds,
                FILE *fout)
{
        int i,j,k,kn;
        int count=0,Mcnt=0;
        float **param;
        float *cent,*expn,*refl,*cont;
        float ttcent,ttrefl,ttexpn,ttcont;
        int maxcount=50;
        float TT[3],T0;
        float thr;
        float Amin=-25.0, Amax=25.0;
        float L,L1;
        float fastest;
        int missed=0;

        /* le rai est-il contraint par le modele ? */
        missed=0;
        for (i=0; i<p->nseg; i++) if (getvelocity(X[i],Y[i],Z[i],mod,1)<0.0) missed++;
        printf("\nMISSED=%d/%d\n",missed,p->nseg/MISSFRAC);
        if (missed>p->nseg/MISSFRAC) printf("!!!!!!!!!!!!!!!!!!!!!!!! missed=%d/%d\n",
                        missed,p->nseg/MISSFRAC);
        else missed=0;

        param=(float**)malloc(3*sizeof(float*));
        for (i=0; i<3; i++) param[i]=(float*)malloc(2*NH*sizeof(float));
        cent=(float*)malloc(2*NH*sizeof(float));
        refl=(float*)malloc(2*NH*sizeof(float));
        expn=(float*)malloc(2*NH*sizeof(float));
        cont=(float*)malloc(2*NH*sizeof(float));

        /* initialisation */
        thr=0.001;
        for (i=0; i<3; i++)
                for (j=0; j<2*NH; j++) param[i][j]=0.0;

        for (j=0; j<2*NH; j++) cent[j]=refl[j]=expn[j]=cont[j]=0.0;

        for (L1=0.,i=0; i<p->nseg; i++) L1+=ds;
        L=L1;
        T0=TT[0]=compT(p->nseg,X,Y,Z,T,ds,&L,param[0],&(param[0][NH]),NH,NULL,0);
        L1=L;
        fprintf(stderr,"L=%10.3f\t",L1);

        for (i=0; i<3; i++) {
                TT[i]=compT(p->nseg,X,Y,Z,T,ds,&L,param[i],&(param[i][NH]),NH,mod,0);
                L=L1;
        }
        T0=TT[0];
        sorttab(3, TT, param);
        fprintf(stderr,"N=%d \n",p->nseg);
        fprintf(stderr,"T=%10.3f\t",T0);
        fprintf(fout,"T= %10.3f\t",T0);

        do {
                int pp=1;
                fastest=TT[0];
                if (missed>p->nseg/MISSFRAC) break;
                for (kn=k=0; k<NH; pp*=-1) {
                        if (pp>0) kn=k;
                        else { kn=k+NH; k++; }

                        memcpy(param[1],param[0],2*NH*sizeof(float));
                        memcpy(param[2],param[0],2*NH*sizeof(float));
                        param[1][kn]=param[0][kn]+Amin;
                        param[2][kn]=param[0][kn]+Amax;
                        count=0;
                        for (i=0; i<3; i++) {
                                TT[i]=compT(p->nseg,X,Y,Z,T,ds,&L,param[i],&(param[i][NH]),NH,mod,0);
                                L=L1;
                        }
                        do {
                                sorttab(3, TT, param);
                                cent[kn]=0.5*(param[0][kn]+param[1][kn]);
                                ttcent=compT(p->nseg,X,Y,Z,T,ds,&L,cent,&(cent[NH]),NH,mod,0); L=L1;

                                refl[kn]=2.0*cent[kn]-param[2][kn];
                                ttrefl=compT(p->nseg,X,Y,Z,T,ds,&L,refl,&(refl[NH]),NH,mod,0); L=L1;

                                if (ttrefl<TT[0]) {
                                        expn[kn]=3.0*cent[kn]-2.0*param[2][kn];
                                        ttexpn=compT(p->nseg,X,Y,Z,T,ds,&L,expn,&(expn[NH]),NH,mod,0); L=L1;
                                        if (ttexpn<TT[0]) {
                                                swtch(&(param[0]),&expn,&(TT[0]),&ttexpn);
                                                swtch(&(param[1]),&expn,&(TT[1]),&ttexpn);
                                                swtch(&(param[2]),&expn,&(TT[2]),&ttexpn);
                                        } else {
                                                swtch(&(param[0]),&refl,&(TT[0]),&ttrefl);
                                                swtch(&(param[1]),&refl,&(TT[1]),&ttrefl);
                                                swtch(&(param[2]),&refl,&(TT[2]),&ttrefl);
                                        }