/* * This file is part of tomo3d * * Copyright (C) 2002, 2003, Sebastien Judenherc <sebastien.judenherc@na.infn.it> * * 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);}