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

#define RE 6370.

int main(int argc, char **argv)
{ struct block_t * blk=NULL;
  int bl,nbb,*hit,maxhit,layer1,nbb0,i,layer,maxlayer=100;
  float X;
  char blkname[255],parname[255],numname[255],
       xyzname[255],hitname[255];
  FILE *f,*fnum,*fxyz,*fhit;
  argv=parseoptions(&argc,argv);
  if (argc<0) { argc*=-1; exit(0); }
  if (argc<3) {
     fprintf(stderr,"%s name maxhit\n",argv[0]);
     exit(0);
  }
  if (!sscanf(argv[2],"%d",&maxhit)) {
     fprintf(stderr,"arg error\n");
     exit(0);
  }
  if (argc>3) {
     if (!sscanf(argv[3],"%d",&maxlayer)) {
        fprintf(stderr,"arg3 error\n");
        exit(0);
     }
  }
  fprintf(stderr,"================== reblk ==================\n");
  sprintf(hitname,"%s.hit",argv[1]);
  sprintf(blkname,"%s.blk",argv[1]);
  sprintf(parname,"%s.par",argv[1]);
  sprintf(numname,"%s.num",argv[1]);
  sprintf(xyzname,"%s.xyz",argv[1]);
  if ((blk=readblock(blkname, &nbb))==NULL) exit(0);
  hit=(int*)calloc(nbb,sizeof(int));
  if ((fnum=fopen(numname,"rt"))==NULL) {
     perror(numname);
     exit(0);
  }
  if ((fxyz=fopen(xyzname,"wt"))==NULL) {
     perror(xyzname);
     exit(0);
  }
  if ((fhit=fopen(hitname,"rt"))==NULL) {
     perror(hitname);
     exit(0);
  }
  fscanf(fnum,"%d",&bl);
  fscanf(fhit,"%d",&(hit[bl]));
  while (!feof(fnum)) {
        fscanf(fnum,"%d",&bl);
        fscanf(fhit,"%d",&(hit[bl]));
  }
  fclose(fnum);
  fclose(fxyz);
  fclose(fhit);
  for (layer1=0,X=blk[0].z; layer1<nbb; layer1++) {
      if (blk[layer1].z>X+blk[0].z/100.0) break;
         X=blk[layer1].z;
      }
  fprintf(stderr,">>>> Layer 1 ending at %d\n",layer1);
  nbb0=nbb;
  layer=0;
  fprintf(stderr,">>>> Entering layer ");
  for (i=layer1; i<nbb0; i++) {
      if (blk[i].z>blk[i-1].z+0.5) {
         layer++;
         fprintf(stderr," %d",layer);
      }
      if (layer>maxlayer) { fprintf(stderr," STOP\n"); break; }
      if (hit[i]>maxhit) {
         blk=(struct block_t*)realloc(blk,(3+nbb)*sizeof(struct block_t));

         memcpy(&(blk[nbb]),&(blk[i]),sizeof(struct block_t));
         blk[nbb].x+=0.5*blk[i].dx;
         blk[nbb].dx=0.5*blk[i].dx;
         blk[nbb].dy=0.5*blk[i].dy;
         nbb++;
         
         memcpy(&(blk[nbb]),&(blk[i]),sizeof(struct block_t));
         blk[nbb].y+=0.5*blk[i].dy;
         blk[nbb].dx=0.5*blk[i].dx;
         blk[nbb].dy=0.5*blk[i].dy;
         nbb++;
         
         memcpy(&(blk[nbb]),&(blk[i]),sizeof(struct block_t));
         blk[nbb].y+=0.5*blk[i].dy;
         blk[nbb].x+=0.5*blk[i].dx;
         blk[nbb].dx=0.5*blk[i].dx;
         blk[nbb].dy=0.5*blk[i].dy;
         nbb++;

         blk[i].dx*=0.5;
         blk[i].dy*=0.5;
      }
  }
  if (nbb==nbb0) {
     fprintf(stderr,"NOTHING TO DO!\n");
     exit(1);
  }
  fprintf(stderr,"Sorting blocks... ");
  qsort((void*)blk, nbb, sizeof(struct block_t),blkcmp);
  f=fopen(blkname,"wt");
  printblocks(f,blk,nbb);
  fclose(f);
  fprintf(stderr,"old: %d new: %d\n",nbb0, nbb);
  return(0);
}