normray.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* NORMRAY: $Test Release: 1.1 $ ; $Date: 1997/10/29 17:23:12 $	*/

#include "par.h"
#include "Triangles/tri.h"
#include "Triangles/sloth.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"                                              				",
" NORMRAY - dynamic ray tracing for normal incidence rays in a sloth model",
"                                              				",
"    normray <modelfile >rayends [optional parameters]			",
"									",
" Optional Parameters:							",
" caustic	 0: show all rays 1: show only caustic rays		",
" nonsurface	 0: show rays which reach surface 1: show all rays      ",
" surface	 0: shot ray from subsurface 1: from surface               ",
" nrays 	 number of location to shoot rays                       ",
" dangle 	 increment of ray angle for one location                ",
" nangle 	 number of rays shot from one location                  ",
" ashift 	 shift first taking off angle                           ",
" xs1 	         x of shooting location                                 ",
" zs1 	         z of shooting location                                 ",
" nangle=101     number of takeoff angles				",
" fangle=-45     first takeoff angle (in degrees)			",
" rayfile        file of ray x,z coordinates of ray-edge intersections	",
" nxz=101        number of (x,z) in optional rayfile (see notes below)	",
" wavefile       file of ray x,z coordinates uniformly sampled in time	",
" nt=101         number of (x,z) in optional wavefile (see notes below)	",
" infofile       ASCII-file to store useful information 		",
" fresnelfile    used if you want to plot the fresnel volumes. 		",
"                default is <fresnelfile.bin> 				",
" outparfile     contains parameters for the plotting software. 	",
"                default is <outpar> 					",
" krecord        if specified, only rays incident at interface with index",
"                krecord are displayed and stored			",
" prim           =1, only single-reflected rays are plotted 		",     
"                =0, only direct hits are displayed  			",
" ffreq=-1       FresnelVolume frequency 				",
" refseq=1,0,0   index of reflector followed by sequence of reflection (1)",
"                transmission(0) or ray stops(-1).			",
"                The default rayend is at the model boundary.		",
"                NOTE:refseq must be defined for each reflector		",
" NOTES:								",
" The rayends file contains ray parameters for the locations at which	",
" the rays terminate.  							",
"									",
" The rayfile is useful for making plots of ray paths.			",
" nxz should be larger than twice the number of triangles intersected	",
" by the rays.								",
"									",
" The wavefile is useful for making plots of wavefronts.		",
" The time sampling interval in the wavefile is tmax/(nt-1),		",
" where tmax is the maximum time for all rays.				",
"									",
" The infofile is useful for collecting information along the		",
" individual rays. The fresnelfile contains data used to plot 		",
" the Fresnel volumes. The outparfile stores information used 		",
" for the plotting software.						",
"									",
NULL};

/*
 * AUTHOR:  Dave Hale, Colorado School of Mines, 02/16/91
 * MODIFIED:  Andreas Rueger, Colorado School of Mines, 08/12/93
 *	Modifications include: functions writeFresnel, checkIfSourceIsOnEdge;
 *		options refseq=, krecord=, prim=, infofile=;
 *		computation of reflection/transmission losses, attenuation.
 * MODIFIED: Boyi Ou, Colorado School of Mines, 4/14/95
 */
/* Notes:
 This code can shoot rays from specified interface by users, normally you
need to use gbmodel2 to generate interface parameters for this code, both
code have a parameter named nrays, it should be same. If you just want to
shoot rays from one specified location, you need to specify xs1,zs1,
otherwise, leave them alone. If you want to shoot rays from surface, you need
to define surface equal to 1. The rays from one location will be
approximately symetric with direction Normal_direction - ashift.(if nangle is
odd, it is symetric, even, almost symetric. The formula for the first take
off angle is: angle=normal_direction-nangle/2*dangle-ashift. If you only want to
see caustics, you specify caustic=1, if you want to see rays which does not
reach surface, you specify nonsurface=1. 
*/
/**************** end self doc ***********************************/


#define TAPER 0.9

typedef struct RSStruct {
	float sigma,x,z,px,pz,t;
	float q1,p1,q2,p2;
	int kmah,nref;
	EdgeUse *eu;
        float atten;
        float ampli,ampliphase;
	Face *f;
	struct RSStruct *rsNext;
} RayStep;

typedef struct RCStruct {
	int k;
	int nhits;
	int *hitseq;
} RefCheck;

/* prototypes for functions defined and used internally */
void shootRays (Model *m, float xs, float zs,
	int nangle, float dangle, float fangle,
	RefCheck *rc, int kk,
	RayStep *rsp[], RayEnd re[], FILE *infofp);
void writeRays (FILE *rayfp, int nxz, int nray, RayStep *rs[], RayEnd  
	re[], int krecord, int prim, FILE *infofp);
void writeFresnel (Model *m,FILE *rayfp,FILE *infofp, int nxz, int nray,
	RayStep *rs[], RayEnd re[], int krecord, float freq, int prim,
	FILE *fresnelfp);
void writeWaves (FILE *wavefp, int nt, int nray, RayStep *rs[]);
static RayStep* traceRayInTri (RayStep *rs);
static RayStep* traceRayAcrossEdge (RayStep *rs,FILE *infofp);
static RayStep* reflectRayFromEdge (RayStep *rs,FILE *infofp);
static void evaluateDynamic (float sigma, float px, float pz, float dsdx,
	float dsdz, float *q1, float *p1, float *q2, float *p2, int *kmah);
static void evaluateDynamic2 (float sigma, float px, float pz, float dsdx,
	float dsdz, float *q1, float *p1, float *q2,float *p2,
	float s0, float s1, float s2); 
int checkIfSourceOnEdge(Face *tris, float zs, float xs);

int caustic;
int raycount; 
int nonsurface;
int surface;
int nrays;

/* the main program */
int main (int argc, char **argv)
{
	int i,ii,ia,nangle,nxz,nt,kk,nseq,krecord,prim;
	int **hitseqa=NULL,*nrefseq=NULL,ibuf=2;
	float *xs,*zs,*as,fangle,dangle,xs1,zs1,ffreq,ashift;
	char *rayfile,*wavefile,*infofile,
		*fresnelfile="fresnelfile.bin",*outparfile="outpar";
	Model *m;
	FILE *rayfp=NULL,*wavefp=NULL,*infofp=NULL;
	FILE *outparfp=NULL,*fresnelfp=NULL;
	FILE *xfp,*zfp,*afp=NULL;

	RayStep **rsp=NULL;
	RayEnd *re=NULL;	
	RefCheck *rc;

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);
	
	/* read model */
	m = readModel(stdin);

	/* get optional parameters */
	
	if (!getparint("nrays",&nrays)) nrays=1;
	if (!getparint("caustic",&caustic)) caustic=0;
	if (!getparint("surface",&surface)) surface=0;
	if (!getparint("nonsurface",&nonsurface)) nonsurface=0;
	if (!getparint("nangle",&nangle)) nangle = 1;
	if (!getparint("nxz",&nxz)) nxz = 101;
	if (!getparint("nt",&nt)) nt = 101;
	if (!getparint("krecord",&krecord)) krecord = INT_MAX;
	if (!getparint("prim",&prim)) prim = INT_MAX;
	if (!getparfloat("ffreq",&ffreq)) ffreq = -1;
	if (!getparfloat("dangle",&dangle)) dangle = 0;
	if (!getparfloat("ashift",&ashift)) ashift = 0;
	if (!getparfloat("xs1",&xs1)) xs1 = 0;
	if (!getparfloat("zs1",&zs1)) zs1 = 0;

	/* open files as necessary */
	if (getparstring("rayfile",&rayfile)) rayfp = efopen(rayfile,"w");
	if (getparstring("wavefile",&wavefile)) wavefp = efopen(wavefile,"w");
	if (getparstring("infofile",&infofile)) infofp = efopen(infofile,"w");
	getparstring("fresnelfile",&fresnelfile);
	if (ffreq>0) fresnelfp = efopen(fresnelfile,"w");
	getparstring("outparfile",&outparfile);
	outparfp = efopen(outparfile,"w");

	/* number of reflection sequences */
        kk = countparname("refseq");

        /* error checks */
	if(prim!=INT_MAX && krecord==INT_MAX)
		err("krecord= must be specified with prim=\n");
	if(prim!=INT_MAX && kk==0)
		err("must define target reflector(s) via refseq= \n");

	/* allocate space for raychecks */
	rc = (RefCheck*)ealloc1(kk,sizeof(RefCheck));

	/* allocate space for hitseq */
	if (kk>0) {
		hitseqa = (int**)ealloc1(kk,sizeof(int*));
		nrefseq = ealloc1int(kk);
	}

	/* determine reflection sequences */
	for (i=0; i<kk; ++i) {
		nseq = countnparval(i+1,"refseq");
		nrefseq[i] = nseq;
		hitseqa[i] = ealloc1int(nrefseq[i]+ibuf);
		getnparint(i+1,"refseq",hitseqa[i]);
		rc[i].k = hitseqa[i][0];
		rc[i].nhits = 0;
		rc[i].hitseq = hitseqa[i]+1;
	}

	/* print reflection sequence information if requested */
	if (infofp!=NULL) {
		fprintf(infofp,
			"THIS FILE CONTAINS RAY TRACING INFORMATION\n\n"
			"*********************************************\n"
			"defined reflection/transmission sequences:\n"
			"*********************************************\n");
		for (i=0; i<kk; ++i) {
			fprintf(infofp,"interface %d:",rc[i].k);
			for(ii=0; ii<nrefseq[i]-1; ++ii)
				fprintf(infofp," %d%s",
					rc[i].hitseq[ii],
					(ii<nrefseq[i]-2)?",":"");
			fprintf(infofp,"\n");
		}
		fprintf(infofp,"Interfaces without defined "
				"refseq are transmitting.\n");
	}
	

	/*allocate space to xs, zs, as */
	xs=alloc1float(nrays);
        zs=alloc1float(nrays);
        as=alloc1float(nrays);

	/* open file of XInterface,ZInterface,AInterface*/
        xfp=fopen("XInterface","r");
        zfp=fopen("ZInterface","r");

        if(surface==0)
         afp=fopen("AInterface","r");

	fread(xs,sizeof(float),nrays,xfp);
        fread(zs,sizeof(float),nrays,zfp);


	if(surface==0)
           fread(as,sizeof(float),nrays,afp);
        else
           for(ia=0;ia<nrays;++ia)
               as[ia]=0;


	/* convert angles to radians and determine angle increment */

        dangle *= PI/180;
        ashift *= PI/180;

        if(nrays==1){
          xs[0]=xs1;
          zs[0]=zs1;
          as[i]=180;
        }

	for(i=0;i<nrays;++i){

	/* allocate space for ray step lists and ray ends */
	rsp = (RayStep**)ealloc1(nangle,sizeof(RayStep*));
	re = (RayEnd*)ealloc1(nangle,sizeof(RayEnd));

	fangle = as[i]*PI/180 - nangle/2*dangle -ashift;

	/* shoot rays */
	shootRays(m,xs[i],zs[i],nangle,dangle,fangle,rc,kk,rsp,re,infofp);
	
	/* if Fresnel Volumes are requested */
	if (fresnelfp!=NULL && rayfp!=NULL) {

		/* write Fresnel Volumes and associated rays */
		writeFresnel(m,rayfp,infofp,nxz,nangle,rsp,re,krecord,
			ffreq,prim,fresnelfp);

	} else if (rayfp!=NULL) {

		/* if requested, write rays to rayfile */
		writeRays(rayfp,nxz,nangle,rsp,re,krecord,prim,
			infofp);
	}

	}


	fprintf(outparfp,"%d\n",raycount);

	/* if requested, write waves to wavefile */
	if (wavefp!=NULL) writeWaves(wavefp,nt,nangle,rsp);

	/* write ray ends */
	if (efwrite(re,sizeof(RayEnd),nangle,stdout)!=nangle)
		err("Error writing ray ends to stdout!\n");

	return EXIT_SUCCESS;
}

void shootRays (Model *m, float xs, float zs,
	int nangle, float dangle, float fangle,
	RefCheck *rc, int kk,
	RayStep *rsp[], RayEnd re[], FILE *infofp)
/*****************************************************************************
Shoot rays via dynamic ray tracing for a sloth model
******************************************************************************
Input:
m		trianglulated sloth model
xs		horizontal coordinate of source
zs		vertical coordinate (depth) of source
nangle		number of ray takeoff angles
dangle		increment in ray takeoff angle
fangle		first ray takeoff angle
rc		pointer to reflector structure
kk		number of reflectors

Output:
rsp		array[nangle] of pointers to linked list of RaySteps
re		array[nangle] of RayEnds
******************************************************************************
Notes:
Rays are traced from triangle to triangle.
Whenever a ray crosses a triangle edge, two RaySteps are added
to the linked list of RaySteps, one for each side of the edge.
The RaySteps are useful if ray parameters along the entire raypath
are required.  The RayEnds are useful if parameters at the ray
endpoint are required.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 02/28/91
Modified: Andreas Rueger, Colorado School of Mines, 08/12/93
******************************************************************************/
{
	int iangle,refl,stop,i,temp,test;
	float angle,sb,zmax,xmax,eps;
	Tri *tris;
	FaceAttributes *fa;
	EdgeAttributes *ea;
	RayStep *rs,*rslast;

	/* determine triangle in model containing source location */
	tris = insideTriInModel(m,NULL,zs,xs);
	
        /* if source sits on edge, move it by eps */
        if(((test=checkIfSourceOnEdge(tris,zs,xs)) != 0)) {
		zmax = m->xmax;
		eps = m->eps;
		if (test==2) eps = 5*eps;
		zs = (zs<zmax-eps) ? zs+eps : zs-eps;

		/* if we are still on an edge, move x value */
		if(checkIfSourceOnEdge(tris,zs,xs) != 0) {
			xmax = m->ymax;
			xs = (xs<xmax-eps) ? xs+eps : xs-eps;
		}

		/* find new triangle */
		tris = insideTriInModel(m,NULL,zs,xs);
	}

	/* determine sloth at source location (beginning of ray) */
	fa = tris->fa;
	sb = fa->s00+fa->dsdx*xs+fa->dsdz*zs;

	/* loop over takeoff angles */
	for (iangle=0,angle=fangle; iangle<nangle; ++iangle,angle+=dangle) {

		/* initialize linked list of ray steps */
		rs = rsp[iangle] = (RayStep*)ealloc1(1,sizeof(RayStep));
		rs->sigma = 0.0;
		rs->x = xs;
		rs->z = zs;
		rs->px = sin(angle)*sqrt(sb);
		rs->pz = cos(angle)*sqrt(sb);
		rs->q1 = 1.0;
		rs->p1 = 0.0;
		rs->q2 = 0.0;
		rs->p2 = 1.0;
		rs->t = 0.0;
		rs->kmah = 0;
		rs->nref = 0;
		rs->eu = NULL;
		rs->f = tris;
		rs->ampli = 1;
		rs->ampliphase = 0;
		rs->atten = 0;
		rs->rsNext = NULL;

		if (infofp!=NULL)
			fprintf(infofp,"\n"
				"****************************************\n"
				"RAY WITH TAKEOFF ANGLE: %g (degrees)\n"
				"****************************************\n",
				angle*180.0/PI);

		/* trace ray */
		do {
			/* trace ray in triangle */
			rs = traceRayInTri(rs);

			/* remember last ray step */
			rslast = rs;

			/* edge attributes */
			ea = rs->eu->e->ea;

			if (ea==NULL) {

				/* null edges are transmitting */
				refl = stop = temp = 0;

			} else {

				/* check if edge is reflecting or stopping */
				for (i=0,temp=ea->k,refl=stop=0; i<kk; ++i) {
					if (temp==rc[i].k) {
						if (rc[i].hitseq[0]==1)
							refl = 1;
						if (rc[i].hitseq[0]==-1)
							stop = 1;
						++(rc[i].nhits);
						++(rc[i].hitseq);
						break;
					}
				}
			}

			if (infofp!=NULL)
				fprintf(infofp,"^^^\n"
					"Interaction with interface %d at "
					"(x=%g,z=%g).\n",
					temp,rslast->x,rslast->z);

			/* if ray at stopping edge, stop */
			if (stop) {
				if (infofp!=NULL)
					fprintf(infofp,"Hit stopping edge.  "
							"Ray stopped.\n");
				break;
			}

			/* else if ray at reflecting edge, reflect */
			else if (refl)
				rs = reflectRayFromEdge(rs,infofp);

			/* else trace ray across edge */
			else
				rs = traceRayAcrossEdge(rs,infofp);
			
		} while (rs!=NULL);
		
		/* reinitialize RayChecks for new ray */
		for (i=0; i<kk; ++i) {
			rc[i].hitseq = rc[i].hitseq-rc[i].nhits;
			rc[i].nhits = 0;
		}

		/* save ray end parameters */
		fa = rslast->f->fa;
		re[iangle].sigma = rslast->sigma;
		re[iangle].x = rslast->x;
		re[iangle].z = rslast->z;
		re[iangle].px = rslast->px;
		re[iangle].pz = rslast->pz;
		re[iangle].t = rslast->t;
		re[iangle].q1 = rslast->q1;
		re[iangle].p1 = rslast->p1;
		re[iangle].q2 = rslast->q2;
		re[iangle].p2 = rslast->p2;
		re[iangle].kmah = rslast->kmah;
		re[iangle].nref = rslast->nref;
		re[iangle].sb = sb;
		re[iangle].se = fa->s00+fa->dsdx*rslast->x+fa->dsdz*rslast->z;
		re[iangle].dsdxe = fa->dsdx;
		re[iangle].dsdze = fa->dsdz;
		re[iangle].ab = angle;
		re[iangle].kend = temp;
		re[iangle].ampli = rslast->ampli;
		re[iangle].ampliphase = rslast->ampliphase;
		re[iangle].atten = rslast->atten;
		re[iangle].dangle = dangle;

		/* optional output of rayend information */
		if (infofp!=NULL) {
			fprintf(infofp,
				"---------------------------------------"
				"------------\n\n"
				"The following information is stored at "
				"the rayend:\n\n");
			fprintf(infofp,
				"Ray stops at (%g,%g) at interface %d\n",
				re[iangle].x,re[iangle].z,re[iangle].kend);
			fprintf(infofp,
				"Takeoff angle=%g   Emergence angle=%g\n",
				re[iangle].ab*180./PI,-asin(re[iangle].px
				*1./sqrt(re[iangle].se))*180./PI);
			fprintf(infofp,"Takeoff velocity=%g   "
				"Emergence velocity=%g\n",
				1./sqrt(sb),1./sqrt(re[iangle].se));
			fprintf(infofp,"Slowness components px=%g   pz=%g\n",
				re[iangle].px,re[iangle].pz);
			fprintf(infofp,"Traveltime=%g   Sigma=%g\n",
				re[iangle].t,re[iangle].sigma);
			fprintf(infofp,"kmah index=%d   "
				"Number of reflections=%d\n",
				re[iangle].kmah,re[iangle].nref);
			fprintf(infofp,"Ray propagator q1=%g   q2=%g\n",
				re[iangle].q1,re[iangle].q2);
			fprintf(infofp,"Ray propagator p1=%g   p2=%g\n",
				re[iangle].p1,re[iangle].p2);
			fprintf(infofp,"Attenuation factor=%g   "
				"Angle increment=%g\n",
				re[iangle].atten,dangle*180./PI);
			fprintf(infofp,"Refl/Transm Amplitude=%g   Phase=%g\n",
				re[iangle].ampli,re[iangle].ampliphase); 
			fprintf(infofp,
				"---------------------------------------"
				"------------\n\n");
		}
	}
}

void writeRays (FILE *rayfp, int nxz, int nray, RayStep *rs[],
	RayEnd re[], int krecord, int prim, FILE *infofp)
/* for each ray, write x,z pairs */
{
	int iray,nrs,irs,nxzw,ns,is,icount;
	float sigmamax=0.0,sigma1,sigma2,x,z,px,pz,dsdx,dsdz,
		ds,xs,zs,dsigma;
	RayStep *rsi;
	FaceAttributes *fa;

	/* initialize counter */
	icount = 0;

	/* loop over rays */
	for (iray=0; iray<nray; ++iray) {

		if(caustic==0)
		   re[iray].kmah=1;

		if(nonsurface==1)
		   re[iray].z=0;

		if(re[iray].kmah!=0 && re[iray].z<0.000001){
		if (krecord==INT_MAX || (krecord==re[iray].kend
			&& (prim==INT_MAX || prim==re[iray].nref))) {

			/* count rays */
			++icount;
			++raycount;

			/* count ray steps while determining maximum sigma */
			for (nrs=0,rsi=rs[iray]; rsi!=NULL;
				++nrs,rsi=rsi->rsNext)
				sigmamax = rsi->sigma;

			/* initialize number of (x,z) written for this ray */
			nxzw = 0;

			/* loop over ray steps */
			for (irs=0,rsi=rs[iray]; irs<nrs;
				++irs,rsi=rsi->rsNext) {

				/* if sufficient number of (x,z) */