sutetraray.c

su 的源代码库

      if (nygd<2) err("nygd too small");
      nygd=(nygd/2)*2+1;

      if (!getparfloat("fygd",&fygd)) fygd=ymin;
      if (fygd<ymin) err("fygd too small");
      if (!getparfloat("dygd",&dygd)) dygd=(ymax-ymin)/(nygd-1);
      if (fygd+(nygd-1)*dygd>ymax*1.01) err("fygd+(nygd-1)*dygd too large");

      if (verbose) {
	    fprintf(jpfp,"shootupward=%d\n",shootupward);
            fprintf(jpfp,"nxgd=%d\nnygd=%d\n",nxgd,nygd);
            fprintf(jpfp,"dxgd=%f\ndygd=%f\n",dxgd,dygd);
            fprintf(jpfp,"fxgd=%f\nfygd=%f\n",fxgd,fygd);
      }

      if (getparstring("rayfile",&rayfile))
	    if ((rayfp=fopen(rayfile,"w"))==NULL)
	      	  err("Can not open rayfile to write");

      if (getparstring("wffile",&wffile))
            if ((wffp=fopen(wffile,"w"))==NULL)
                  err("Can not open wffile to write");

      if (wffp!=NULL) {
	    if (!getparint("ifwf2dump",&ifwf2dump)) ifwf2dump=20;
            if (!getparint("nwf2dump",&nwf2dump)) nwf2dump=1;
      }

      if (verbose) 
	    fprintf(jpfp,"ifwf2dump=%d, nwf2dump=%d\n",
		  ifwf2dump,nwf2dump);

      if (getparstring("crfile",&crfile))
            if ((crfp=fopen(crfile,"w"))==NULL)
                  err("Can not open crfile to write");

      if (getparstring("sttfile",&sttfile))
            if ((sttfp=fopen(sttfile,"w"))==NULL)
                  err("Can not open sttfile to write");

      if (!getparint("maxntris",&maxntris))
	    maxntris=1000;
	
      if (!getparint("maxnrays",&maxnrays))
            maxnrays=600;

      if (!getparint("maxnsegs",&maxnsegs))
            maxnsegs=100;

      if (!getparfloat("edgemax",&edgemax))
	    edgemax=4;

      if (!getparfloat("fsx",&fsx)) fsx=fxgd+nxgd/2.0*dxgd;
      if (!getparfloat("fsy",&fsy)) fsy=fygd+nygd/2.0*dygd;

      if (!getparint("nsx",&nsx)) nsx=1;
      if (!getparint("nsy",&nsy)) nsy=1;
      if (!getparint("nsz",&nsz)) nsz=1;

      ns=nsx*nsy;

      if (!getparint("nxt",&nxt)) nxt=(nxgd/2)*2+1;
      if (!getparint("nyt",&nyt)) nyt=(nygd/2)*2+1;
      nxt=(nxt/2)*2+1;
      nyt=(nyt/2)*2+1;

      if (!getparfloat("dsx",&dsx)) dsx=dxgd;
      if (!getparfloat("dsy",&dsy)) dsy=dygd;
      if (!getparfloat("dsz",&dsz)) dsz=(zmax-zmin)/2;

      ixsf=(fsx-fxgd)/dxgd;
      if (ixsf<nxt/2) err("ixsf=%d should >= nxt/2=%d. Use larger fsx or smaller nxt",
	    ixsf,nxt/2);

      ixsr=MAX(1,dsx/dxgd);
      if (ixsf+ixsr*nsx+nxt/2 > nxgd) 
            err("ixsf+ixsr*nsx+nxt/2=%d should <= nxgd=%d",
                  ixsf+ixsr*nsx+nxt/2,nxgd);

      iysf=(fsy-fygd)/dygd;
      if (iysf<nyt/2) err("iysf=%d should >= nyt/2=%d. Use larger fsy or smaller nyt",
            iysf,nyt/2);

      iysr=MAX(1,dsy/dygd);
      if (iysf+iysr*nsy+nyt/2 > nygd) 
            err("iysf+iysr*nsy+nyt/2=%d should <= nygd=%d",
                  iysf+iysr*nsy+nyt/2,nygd);

      fsx=fxgd+ixsf*dxgd;
      fsy=fygd+iysf*dygd;
      fsx+=0.01*dxgd;
      fsy+=0.01*dygd;

      dsx=ixsr*dxgd;
      dsy=iysr*dygd;

      if (!getparfloat("fsz",&fsz)) fsz=zmin+(zmax-zmin)/2.0;

      if (!getparint("ntakeoff",&ntakeoff)) ntakeoff=5;
      if (!getparint("nazimuth",&nazimuth)) nazimuth=15;

      if (!getparfloat("tmax",&tmax)) tmax=10.0;

      if (!getparfloat("dtwf",&dtwf)) dtwf=0.05;

      if (!getparint("ntwf",&ntwf)) {
            ntwf=(int)(tmax/dtwf)+1;
            ntwf=MAX(200,ntwf);
      }

      if (ntwf<ifwf2dump+nwf2dump) {
            fprintf(jpfp,"ntwf < ifwf2dump+nwf2dump, reset=%d\n",
                  ifwf2dump+nwf2dump);
            ntwf=ifwf2dump+nwf2dump;
      }

      if (ntwf>2000) err("ntwf>2000, too large");

      if (verbose) {
	    fprintf(jpfp,"nxgd=%d\nnygd=%d\ndxgd=%g\n",
		  nxgd,nygd,dxgd);
	    fprintf(jpfp,"dygd=%g\nntakeoff=%d\nnazimuth=%d\n",
		  dygd,ntakeoff,nazimuth);
            fprintf(jpfp,"takeoff=%f\n",takeoff);
	    fprintf(jpfp,
		  "fsx=%g,dsx=%g,nsx=%d\nfsy=%g,dsy=%g,nsy=%d\nfsz=%g,dsz=%g,nsz=%d\n",
		  fsx,dsx,nsx,fsy,dsy,nsy,fsz,dsz,nsz);
	    fprintf(jpfp,"nxt=%d\nnyt=%d\n",nxt,nyt);
	    fprintf(jpfp,"maxntris=%d\nmaxnrays=%d\nmaxnsegs=%d\nedgemax=%g\n",
		  maxntris,maxnrays,maxnsegs,edgemax);
	    fprintf(jpfp,"ntwf=%d,dtwf=%g,tmax=%g\n",ntwf,dtwf,tmax);
      }

      if (nxt<0 || nxt>nxgd) err("nxt must be >=0 <=nxgd");
      if (nyt<0 || nyt>nygd) err("nyt must be >=0 <=nygd");

      if ((fsx+(nsx-1)*dsx>fxgd+(nxgd-1)*dxgd) ||
          (fsy+(nsy-1)*dsy>fygd+(nygd-1)*dygd) ||
	  (fsx<fxgd) || (fsy<fygd)) 
	    err("Source out of medium");

      /* number of sources */
      tro.nhs=nsx;
      tro.nvs=nsy;
      tro.duse=nsz;

      /* first source */
      tro.sdel=ixsf;
      tro.gdel=iysf;
      tro.sdepth=(int)(fsz*1000.0);

      /* source spacing */
      tro.swdep=ixsr;
      tro.gwdep=iysr;
      tro.ep=(int)(dsz*1000.0);

      /* first sample */
      tro.f1=fxgd;
      tro.f2=fygd;

      /* model size: nxgd x nygd */
      tro.gelev=nygd;
      tro.selev=nxgd;

      /* traveltime sample spacing */
      tro.d1=dxgd;
      tro.d2=dygd;

      /* number of traces and samples in each trace */
      tro.ns=nxt;
      tro.ntr=nyt;

      fscanf(stdin,"%d=npoint : x[3] : s : n[3]\n",&npoint);
      fprintf(jpfp,"npoint=%d\n",npoint);
      point=(struct POINT *)malloc(npoint*sizeof(struct POINT));   
      read_point(
      point,
      npoint);
      fprintf(jpfp,"point %d read\n",npoint);

      fscanf(stdin,"%d=nfacet : area : ct[10] : ip[3] : cn[3] : itetra[2]\n",
	    &nfacet);
      fprintf(jpfp,"nfacet=%d\n",nfacet);
      facet=(struct FACET *)malloc(nfacet*sizeof(struct FACET));      
      read_facet(
            facet,
            nfacet);
            fprintf(jpfp,"facet %d read\n",nfacet);

      fscanf(stdin,"%d=ntetra : v : ip[4] : ifacet[4] : ireg : gs[3]\n",&ntetra);
      fprintf(jpfp,"ntetra=%d\n",ntetra);
      tetra=(struct TETRA *)malloc(ntetra*sizeof(struct TETRA));
      
      read_tetra(
            tetra,
            ntetra);
            fprintf(jpfp,"tetra %d read\n",ntetra);

      if (rayfp!=NULL)
            fprintf(rayfp,"%d =Number of shots\n",ns);

      for (isy=0;isy<nsy;isy++) {
	    source[1]=fsy+isy*dsy;
            iytf=iysf+iysr*isy-nyt/2;
	    for (isx=0;isx<nsx;isx++) {
		  source[0]=fsx+isx*dsx;
                  ixtf=ixsf+ixsr*isx-nxt/2;
                  for (isz=0;isz<nsz;isz++) {
	                source[2]=fsz+isz*dsz;
			if (isz<nsz-1) rayfp0=NULL;
			else rayfp0=rayfp;
			/*******************************************************
                        Usually it is not every good to put the source right on
                        a facet, this is because the raytracer may not stably
                        find the tetrahetron that the source lies in. So addition 
			of EPS will help.
			********************************************************/

			/* shooting one shot */
			oneshot(
			      shootupward,
                              maxnsegs,
                              npoint,
                              point,
                              nfacet,
			      facet,
			      ntetra,
			      tetra,
                              irefseq,/* reflector sequence */
                              nrefseq,/* number of reflectors */
			      jpfp,
			      takeoff,/* max takeoff angle */
			      dtwf,	/* time step */
			      ntwf,	/* maximum number of steps */
		              maxntris,/* maximum number of triangles allowed */
        		      maxnrays,/* maximum number of rays allowed */
        		      edgemax,/* maximum edge length in dxgd */
			      nxt,	/* # of x-samples in ttable for sukdmig3d */
			      nyt,	/* # of y-samples in ttable for sukdmig3d */
			      ixtf,	/* # of samples in ttable on left of source */
			      iytf,	/* # of samples in ttable in front of source */
			      source,   /* source position */
			      ntakeoff, /* number of samples in azimuth */
			      nazimuth, /* number of take-off angles */
			      fxgd,     /* first sample in x */
			      dxgd,	/* sample interval in x */
			      nxgd,	/* number of sampels in x */
			      fygd,	/* first sample in y */
			      dygd,	/* sample interval in y */
			      nygd,	/* number of sampels in y */
			      xmin,ymin,xmax,ymax,zmax,
			      rayfp0,	/* file pointer to rayfile */
			      wffp,
			      crfp,
			      sttfp,
			      ifwf2dump,
                              nwf2dump);

                        fprintf(stderr,"shot (%d %d %d) is written\n",
                              isx,isy,isz);
			if (wffp!=NULL) 
				fclose(wffp);
			wffp=NULL;
		  }
       	    }
      }

      if (rayfp!=NULL) fclose(rayfp);
      if (jpfp!=stderr) fclose(jpfp);

      return EXIT_SUCCESS;
}
void WF_normalize(float *n1,float *n2,float *n3)
/**********************************************************************
WF_normalize - normalize vector n1,n2,n3
**********************************************************************
Inputs/outputs:
float *n1	first component of the vector to be normalized to 1
float *n2       second component of the vector to be normalized to 1
float *n3       the third component of the vector to be normalized to 1
******************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
******************************************************************/
{
        float norm;
        norm=*n1**n1+*n2**n2+*n3**n3;
        norm=sqrt(norm);
        *n1=*n1/norm;
        *n2=*n2/norm;
        *n3=*n3/norm;
}

int    
oneshot(int shootupward, int maxnsegs, int npoint, struct POINT *point,
        int nfacet, struct FACET *facet, int ntetra, struct TETRA *tetra, 
        int *irefseq, int nrefseq, FILE *jpfp, float takeoff, float dtwf,
        int ntwf, int maxntris, int maxnrays, float edgemax, 
        int nxt, int nyt, int ixtf, int iytf, float source[3], int ntakeoff,
        int nazimuth, float fxgd, float dxgd, int nxgd, float fygd, 
        float dygd, int nygd, float xmin, float ymin, float xmax, float ymax,
        float zmax, FILE *rayfp, FILE *wffp, FILE *crfp, FILE *sttfp,
        int ifwf2dump, int nwf2dump)  
/*******************************************************************
Inputs:
int shootupward 1 upward, -1 downward
int maxnsegs    max number of ray segments
int npoint      number of point
struct *point   control points
int nfacet      number of facet
struct *facet   facet information
int ntetra      number of tetra
struct *tetra   tetra information
int *irefseq    reflector sequence
int nrefseq     number of reflectors
FILE *jpfp	file pointer to the job print
float takeoff   maximum takeoff angle
float dtwf        time step for wavefront updating
int ntwf          total number of time steps allowed
int maxntris    maximum number of triangles allowed
int maxnrays    maximum number of rays allowed
float edgemax   maximum edge length in dxgd
int nxt         number of x-samples in ttable for sukdmig3d
int nyt         number of y-samples in ttable for sukdmig3d
int ixtf        number of samples in ttable on left of source
int iytf        number of samples in ttable in front of source
float source[3] source position
int ntakeoff    number of samples in azimuth
int nazimuth    number of take-off angles
float fxgd      first sample in x
float dxgd      sample interval in x
int nxgd        number of sampels in x
float fygd      first sample in y
float dygd      sample interval in y
int nygd        number of sampels in y
FILE *rayfp     file pointer to rayfile
FILE *wffp      file pointer to wavefront file
int ifwf2dump   the first wavefront to dump
int nwf2dump    the number of wavefronts to dump

Output:         none
******************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
******************************************************************/
{
      struct WF *wf;

      /*******************************************************
      The traveltime table start from ixgdl, end with ixgdl+nxt-1 in x 
      and start with iygdf and end with iygdf+nyt-1
      *******************************************************/
      float azimuth;   	/* azimuth angle */
      float sintakeoff;	/* sin(takeoff) */
      float cosazimuth;	/* cos(azimuth) */
      float costakeoff;	/* cos(takeoff) */
      float sinazimuth;	/* sin(azimuth) */
      float s;	        /* slowness */
      int itakeoff;	/* index for azimuth */
      int iazimuth;	/* index for takeoff */

      float dazimuth;	/* azimuth increment */
      float dtakeoff;   /* takeoff increment */
      int ithis;	/* index for this ray */
      int inext;
      int itwf;		/* index for t */
      float distmax=0;	/* the maximum edge of the triangle */
      int iedgemax;	/* the index for the largest triangle */
      int iedgemax0;
      int iedge;	/* index for iedgemax */

      enum ENOUGH_NTWF enough_ntwf=WF_TRUE;

      int i1,i2,i3;	
      int ii1,ii2,ii3;
      int iray;	/* index for rays */