ray2bmatrix.c

射线追踪程序

        ptime(0,1);

#ifdef NTHREADS
        fprintf(stderr,"(NTHREADS=%d) ",NTHREADS);
        threads = (pthread_t *) calloc ( NTHREADS , sizeof (pthread_t ));
        assert(threads);
        pt = (struct pt_param *) calloc ( NTHREADS , sizeof (struct pt_param ));
        assert(pt);
        nt=nbt=0;
#endif
        iclock=time(NULL);
        for (i=0; i<Nray; i++) {
                /* read the ray header */
                fread(&(path[i]),sizeof(struct path_t),1,f);
                p=&(path[i]);
                if (p->nseg>NSEGMAX) {
#ifdef NTHREADS
                        int ii;
                        NSEGMAX=p->nseg+50;
                        for (ii=0; ii<NTHREADS; ii++) {
                                pt[ii].X=realloc(pt[ii].X,sizeof(float)*NSEGMAX);
                                pt[ii].Y=realloc(pt[ii].Y,sizeof(float)*NSEGMAX);
                                pt[ii].Z=realloc(pt[ii].Z,sizeof(float)*NSEGMAX);
                                pt[ii].T=realloc(pt[ii].T,sizeof(float)*NSEGMAX);
                        }
#else
                        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);
#endif
                }

                /* read the ray path and time along segments DS */
#ifdef NTHREADS
                fread(pt[nt].X,sizeof(float),p->nseg,f);
                fread(pt[nt].Y,sizeof(float),p->nseg,f);
                fread(pt[nt].Z,sizeof(float),p->nseg,f);
                fread(pt[nt].T,sizeof(float),p->nseg,f);
#else
                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);
#endif

                /* 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;
                /* the ray is going back from the station towards the source */

#ifdef NTHREADS
                pt[nt].id=i;
                pt[nt].bltmp=bltmp;
                pt[nt].blk=blk;
                pt[nt].nbb=nbb;
                pt[nt].i=i;
                pt[nt].maxdepth=maxdepth;
                pt[nt].p=p;
                assert(pthread_create(&(threads[nt]),NULL,
                                        (void *)thr_processray,(void *) &(pt[nt]))==0);
                nt++;
                nbt++;
                if ((nbt==NTHREADS)||(i+1==Nray)) {
                        for (nt=0; nt<nbt; nt++) {
                                assert(pthread_join(threads[nt],(void *) &status)==0);
                                assert(*status==pt[nt].id);
                                pt[nt].id=-1;
                        }
                        nbt=nt=0;
                }
                if (!(i%1000)) {
                    __lclock=0;
                    PROGRESS(((Nray-i)*(time(NULL)-iclock))/(1+i),
                                    Nray*(time(NULL)-iclock)/(1+i));
                }
#else
                processray(bltmp,blk,nbb,X,Y,Z,T,p,i,maxdepth);
                PROGRESS(((Nray-i)*(time(NULL)-iclock))/(1+i),
                                Nray*(time(NULL)-iclock)/(1+i));
#endif
        }
        ptime(0,2);
        fprintf(stderr,"\nCounting blocks... ");
        for (i=0; i<Nray; i++)
                for (bltmp=&(blt[i]); bltmp; bltmp=bltmp->next)
                        hit[bltmp->num]++;

        fprintf(stderr,"done\n");

        fprintf(stderr,"Reordering blocks...");
        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,"done: %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,"G[i][j]... ");
        fprintf(stderr,"init(G)... ");
        G=newmatrix(Nray,Nblk);
        if (!G) {
                fprintf(stderr,"alloc error\n");
                exit(1);
        }
        fprintf(stderr," computing...");

        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] ",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;
                float x;
                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;
                        x=bltmp->T/d[i];
                        if (ACH) x-=(sG[leve][nb]/(float)nG[leve]);
                        sparseputval(G,i,nb,x);
                }

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

        ptime(0,1);
        fprintf(stderr," done\n");
        if (NULL==(f=fopen(sysname,"wb"))) {
                perror(sysname);
                exit(1);
        }
        writematrix(f,G);
        fprintf(stderr,"Writing [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,"[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);
}