statmat.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 <time.h>
#include <unistd.h>
#include <sismoutil.h>
#include "tomo3d.h"

int main(int argc, char **argv)
{ struct block_t * blk=NULL;
  FILE *f;
  float **G;
  float *d;
  int N,M,nbb;
  int *num,i,j,layer;
  char sysname[255],numname[255],blkname[255];
  double minv,maxv,avgv,stdv,colv;
  argv=parseoptions(&argc,argv);
  if (argc<0) { argc*=-1; exit(0); }
  if (argc<2) {
     fprintf(stderr,"%s name event-file\n",argv[0]);
     exit(0);
  }
  sprintf(sysname,"%s.sys",argv[1]);
  sprintf(numname,"%s.num",argv[1]);
  sprintf(blkname,"%s.blk",argv[1]);
  
  if (NULL==(f=fopen(sysname,"rb"))) {
     perror(sysname);
     exit(1);
  }
  fread(&N,sizeof(int),1,f);
  fread(&M,sizeof(int),1,f);
  d=(float*)calloc(N,sizeof(float));
  G=(float**)calloc(N,sizeof(float*));
  for (i=0; i<N; i++)
      if ((G[i]=(float*)calloc(M,sizeof(float)))==NULL) {
         fprintf(stderr,"malloc error (1:%d)\n",i);
         exit(0);
  }
  for (i=0; i<N; i++)
      if (fread(G[i],sizeof(float),M,f)!=M) {
         fprintf(stderr,"Error (1) reading %s\n",sysname);
         exit(0);
      }
  if (fread(d,sizeof(float),N,f)!=N) {
     fprintf(stderr,"Error (2) reading %s\n",sysname);
     exit(0);
  }
  fclose(f);
  /*
  f=fopen("sys.bin","wb");
  for (i=0; i<N; i++)
      for (j=0; j<M; j++) {
          float titi;
          titi=(float)i;
          fwrite(&titi,sizeof(float),1,f);
          titi=(float)j;
          fwrite(&titi,sizeof(float),1,f);
          fwrite(&(G[i][j]),sizeof(float),1,f);
      }
  fclose(f);
  */
  if (NULL==(f=fopen(numname,"rt"))) {
     perror(numname);
     exit(1);
  }
  if ((blk=readblock(blkname, &nbb))==NULL) exit(0);
  num=(int*)calloc(M,sizeof(int));
  for (i=0; i<M; i++)
      if (!fscanf(f,"%d",&(num[i]))) {
         fprintf(stderr,"Error reading %s\n",numname);
         exit(1);
      }
  fclose(f);
  fprintf(stderr,"Blocks: total=%d, used=%d\n",nbb,M);
  for (layer=0,i=0; i<M; i++) {
      printf("Layer %3d: z=%5.0f..",layer,blk[num[i]].z);
      for (j=i,++i; i<M; i++)
          if (blk[num[i]].z>blk[num[i-1]].z) break;
      printf("%5.0f     N=%3d (%3d,%3d)\n",
             blk[num[j]].z+blk[num[j]].dz,i-j,j,i-1);
      i--;
      layer++;
  }
  minv=maxv=G[0][0];
  avgv=0.0;
  for (i=0; i<N; i++) {
      for (j=0; j<M; j++) {
          if (G[i][j]>maxv) maxv=G[i][j];
          if (G[i][j]<minv) minv=G[i][j];
          avgv+=G[i][j];
      }
  }
  avgv/=(float)(N*M);
  printf("G:         max=%.4e\n           avg=%.4e\n",
         maxv,avgv);
  fflush(stdout);
  stdv=0.;
  for (i=0; i<N; i++) {
      for (j=0; j<M; j++) {
          stdv+=(G[i][j]-avgv)*(G[i][j]-avgv);
      }
  }
  stdv/=(float)(N*M);
  printf("           std=%.4e\n",sqrt(stdv));
  avgv=0.;
  for (j=0; j<M; j++) {
      colv=0.;
      for (i=0; i<N; i++) colv+=G[i][j];
      if (!j) minv=maxv=colv;
      if (minv>colv) minv=colv;
      if (maxv<colv) maxv=colv;
      avgv+=colv;
  }
  avgv/=(float)(M);
  printf("COL(G):    min=%.4e\n           max=%.4e\n           avg=%.4e\n",
         minv,maxv,avgv);
  stdv=0.;
  for (j=0; j<M; j++) {
      colv=0.;
      for (i=0; i<N; i++) colv+=G[i][j];
      stdv+=(colv-avgv)*(colv-avgv);
  }
  stdv/=(float)(M);
  printf("           std=%.4e\n",sqrt(stdv));
  avgv=0.0;
  for (i=0; i<M; i++) {
      colv=0.;
      for (j=0; j<M; j++) colv+=G[j][i]*G[j][i];
      if (!i) maxv=minv=colv;
      if (colv>maxv) maxv=colv;
      if (colv<minv) minv=colv;
      avgv+=colv;
  }
  avgv/=(float)M;
  printf("DIAG(GtG): min=%.4e\n           max=%.4e\n           avg=%.4e\n",
         minv,maxv,avgv);
  fflush(stdout);
  stdv=0.;
  for (i=0; i<M; i++) {
      colv=0.;
      for (j=0; j<M; j++) colv+=G[j][i]*G[j][i];
      stdv+=(colv-avgv)*(colv-avgv);
  }
  stdv/=(float)M;
  printf("           std=%.4e\n",sqrt(stdv));
  f=fopen("block","wt");
  if (!f) {
          perror("block");
          exit(1);
  }
  for (layer=0,i=0; i<M; i++) {
      for (j=i,++i; i<M; i++) {
              if (blk[num[i]].z>blk[num[i-1]].z) break;
      }
      fprintf(f,"%d %d %d\n",layer+1,i-j,i);
      layer++;
  }
  fclose(f);
  return(0);
}