samplegeo.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 OK fprintf(stderr,"FILE %s, LINE %d\n",__FILE__,__LINE__);
#define RE 6370.

struct pt_t { float X,Y,Z;
              float pert,res,sd,hit;
              int lay;
};

int main(int argc, char **argv)
{ struct pt_t *pt=NULL;
  int lay=0,i,Nblk=0,j;
  double pert,sd,res,sum,dhit;
  char slay[255];
  float Z;
  double minlon,maxlon,minlat,maxlat;
  double lat,lon,X,Y;
  float wfact=-1.0,width=-1.0,inc=-1.0;
  struct lambert proj;
  char *basename=NULL;
  char geoname[255],layname[255];
  FILE *flay,*fgeo;
  argv=parseoptions(&argc,argv);
  if (argc<0) { argc*=-1; exit(0); }
  for (i=1; i<argc; i++) {
      if (argv[i][0]=='-') {
         switch (argv[i][1]) {
                case 'n' : if (argc<=i+1) goto error;
                           basename=argv[++i];
                           break;
                case 'I' : if (i+1>=argc) goto error;
                           if (!sscanf(argv[++i],"%f",&inc)) goto error;
                           break;
                case 'W' : if (i+1>=argc) goto error;
                           if (!sscanf(argv[++i],"%f",&wfact)) goto error;
                           break;
                default:   goto usage;
         }
      }
  }
  if (!basename) goto usage;
  if (inc<0.) goto usage;
  if (wfact<0.) goto usage;
  sprintf(geoname,"%s.geo",basename);
  sprintf(layname,"%s.lay",basename);
  flay=fopen(layname,"rt");
  if (!flay) { perror(layname); exit(0); }
  { double lat=48.0,lon=-2.5;
    if (fscanf(flay,"%lf %lf",&lat,&lon)!=2) {
       fprintf(stderr,"ERR: %s\n",layname);
       exit(0);
    }
    initproj(13,lon,lat,0.0,&proj);
  }
  if ((fgeo=fopen(geoname,"rt"))==NULL) {
     perror(geoname);
     exit(0);
  }
  lay=0;
  fscanf(fgeo,"%lf %lf %f %lf %lf %s %lf %lf\n",
               &lon,&lat,&Z,&pert,&sd,slay,&dhit,&res);
  while (!feof(fgeo)) {
        pt=(struct pt_t*)realloc(pt,(Nblk+1)*sizeof(struct pt_t));
        sscanf(slay+5,"%d",&lay);
        geo2lam(lat,lon,&X,&Y,&proj);
        pt[Nblk].pert=pert;
        pt[Nblk].sd=sd;
        pt[Nblk].lay=lay;
        pt[Nblk].hit=dhit;
        pt[Nblk].res=res;
        pt[Nblk].X=X/1000.;
        pt[Nblk].Y=Y/1000.;
        pt[Nblk].Z=Z;
        if (!Nblk) {
           minlat=maxlat=lat;
           minlon=maxlon=lon;
        }
        if (lat<minlat) minlat=lat;
        if (lon<minlon) minlon=lon;
        if (lat>maxlat) maxlat=lat;
        if (lon>maxlon) maxlon=lon;
        fscanf(fgeo,"%lf %lf %f %lf %lf %s %lf %lf\n",
               &lon,&lat,&Z,&pert,&sd,slay,&dhit,&res);
        Nblk++;
  }
  fclose(fgeo);
  fprintf(stderr,"lat: %f -> %f\n",minlat,maxlat);
  fprintf(stderr,"lon: %f -> %f\n",minlon,maxlon);
  fprintf(stderr,"inc=%f\n",inc);
  fprintf(stderr,"Nblk=%d\n",Nblk);

  sprintf(geoname,"%s-s.geo",basename);
  fprintf(stderr,"Output: %s\n",geoname);
  if ((fgeo=fopen(geoname,"wt"))==NULL) {
     perror(geoname);
     exit(0);
  }
  lay++;

  minlat=46.0;
  maxlat=49.5;
  minlon=-4.5;
  maxlon=0.5;
  fscanf(flay,"%*d %lf %*d %lf %*f %*f %*f",&X,&Y);
  fprintf(stderr,"Layer: ");
  for (i=1; i<lay; i++) {
      double dist2;
      int startlay=-1,stoplay=-1;
      for (j=0; j<Nblk; j++) {
          if (i==pt[j].lay)
             if (startlay<0) startlay=j;
          if (i!=pt[j].lay)
             if (startlay>=0) { stoplay=j; break; }
      }
      if ((stoplay<0)&&(startlay>=0)) stoplay=Nblk;
      fscanf(flay,"%*d %lf %*d %lf %*f %*f %*f",&X,&Y);
      width=0.5*wfact*(X+Y);
      fprintf(stderr,"%d(%3d,%3d) [%3.1fkm]\n       ",
              i,startlay,stoplay,width);
      for (lat=minlat; lat<maxlat; lat+=inc)
          for (lon=minlon; lon<maxlon; lon+=inc) {
              geo2lam(lat,lon,&X,&Y,&proj);
              X/=1000.;
              Y/=1000.;
              dhit=sum=sd=res=pert=0.;
              for (j=startlay; j<stoplay; j++) {
                  dist2=(X-pt[j].X)*(X-pt[j].X)+(Y-pt[j].Y)*(Y-pt[j].Y);
                  Z=pt[j].Z;
                  if (dist2>10.*width*width) continue;
                  dist2=exp(-2.0*dist2/(width*width));
                  sum+=dist2;
                  pert+=dist2*pt[j].pert;
                  sd  +=dist2*pt[j].sd  ;
                  dhit +=dist2*pt[j].hit ;
                  res +=dist2*pt[j].res ;
              }
              if (dhit<1.0) continue;
              dhit/=sum;
              sd/=sum;
              pert/=sum;
              res/=sum;
              fprintf(fgeo,"%10e %10e %10e %10e %10e layer%d %10f %10f\n",
                           lon,lat,Z,pert,sd,i,dhit,res);
          }
  }
  fprintf(stderr,"\n");
  fclose(fgeo);
  fclose(flay);
  return(0);
error:
  fprintf(stderr,"%s: bad argument: '%s'\n",argv[0],argv[i]);
usage:
  fprintf(stderr,"%s -n name -I inc -W width\n",argv[0]);
  return(1);
}