1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202

/*
* $Id: track.h,v 1.1 20050918 22:04:56 dhmunro Exp $
* Routines for tracking a straight ray in 3D through a cylindrical mesh.
*/
/* Copyright (c) 2005, The Regents of the University of California.
* All rights reserved.
* This file is part of yorick (http://yorick.sourceforge.net).
* Read the accompanying LICENSE file for details.
*/
#ifndef TRACK_H
#define TRACK_H
#include "bound.h"
/*  */
/* The result of a ray tracking operation is a RayPath: */
typedef struct RayPath RayPath;
struct RayPath {
long maxcuts, ncuts; /* physical lengths, current number zones cut */
long *zone; /* [ncuts] array of zone entered, if last exit
* (into region 0) then zone=0 */
double *ds; /* [ncuts] path length in zone (0.0 if zone=0) */
long *pt1, *pt2; /* [ncuts] arrays, pt1>pt2 is edge of zone cut */
double *f; /* [ncuts] cut at pt1 + (f+0.5)*(pt2pt1) */
double fi,ff; /* fraction of distance into initial (final)
* zone at which tracking started (stopped)
* ds in initial (final) zone is 1fi (1ff) of
* the full length of the ray within the zone */
};
/* Although the mesh lives in a 2D cylindrical coordinate space, a Ray
* is a directed line in Cartesian 3D space. The zaxis is common to
* the two coordinate systems. The xaxis coincides with the raxis
* in a standard (z,r) mesh plot (although, of course, r>=0, while x
* can have either sign). The yaxis comes out of the page of the
* standard mesh plot, toward your face. Hence, (x,y,z) is a right
* handed coordinate frame.
*
* A ray is assumed to lie in a plane of constant y. The angle theta
* is the angle from the +z direction to the ray direction, with
* theta=+pi/2 corresponding to the +xdirection. (This is opposite
* the convention in TDG/DIRT, but corresponds to phi=0 azimuth in
* the usual spherical coordinate system, e.g. Jackson's E&M book.)
*
* The redundant data in the Ray structure allows the ray tracker to
* step from edge to edge without having to compute a lot of extra
* square roots, e.g. to get r=sqrt(x^2+y^2) or sin=sqrt(1cos^2).
* (The alternative approach, used in TDG/DIRT, is to describe the ray
* in an invariant fashion, i.e. independent of the current zone in
* the tracking process.)
*/
typedef struct Ray Ray;
struct Ray {
double cos, sin; /* cos(theta), sin(theta) */
double y, z, x; /* any point on the ray...
* y remains fixed,
* while z and x step from edge to edge */
double r; /* (z,r) are cylindrical coordinates corresponding
* to (x,y,z)  r^2 = x^2 + y^2 */
};
/* Structure to store results of a rayedge intersection calculation
* made by ExitEdge. There are potentially two solutions, the "exit"
* solution, representing a lefttoright crossing in (z,r) space, and
* the "entry" solution, representing a righttoleft crossing.
* Neither solution necessarily exists. Several intermediate results
* of the calculation are saved here, as well as the two roots.
*/
typedef struct RayEdgeInfo RayEdgeInfo;
struct RayEdgeInfo {
double dz, dr; /* edge vector in (z,r) plane */
double area; /* a useful cross product (see EdgeExit) */
double A, B, C; /* Af^2+Bf+C=0 is equation for fx, fn 
* neither B nor C valid unless D>0 */
double D; /* sqrt(discriminant) if >0, else discriminant */
double fx; /* the lefttoright crossing (edge fraction) */
int validx; /* fx not valid unless this is true */
double fn; /* the righttoleft crossing (edge fraction) */
int validn; /* fn not valid unless this is true */
};
/* A ray enters the mesh to begin tracking at one or more EntryPoint's: */
typedef struct EntryPoint EntryPoint;
struct EntryPoint {
EntryPoint *next;
Ray ray;
RayEdgeInfo info;
long zone;
int side;
double f;
double s0;
};
/*  */
/* Given the mesh boundary (i.e. list of boundary edges) and a ray,
* return a linked list of entry point(s), ordered by increasing time.
*/
extern EntryPoint *FindEntryPoints(Boundary *boundary, Ray *ray);
extern void FreeEntryPoints(EntryPoint *entry);
/* Starting at the given entry point(s), track the ray through the mesh.
* The path arrays will be lengthened if necessary, but never shortened.
* sLimits[0] <= s[1]<=s[ncuts2] <= sLimits[1], if sLimits[0]<sLimits[1]
* on entry. (Note that s[0] and s[ncuts1] may lie outside range.)
*/
extern void RayTrack(Mesh *mesh, EntryPoint *entry, RayPath *path,
double *sLimits);
/* Track the ray through the mesh, assuming it represents a sphere.
* The path arrays will be lengthened if necessary, but never shortened.
* sLimits[0] <= s[1]<=s[ncuts2] <= sLimits[1], if sLimits[0]<sLimits[1]
* on entry. (Note that s[0] and s[ncuts1] may lie outside range.)
*/
extern void RayTrackS(Mesh *mesh, Ray *ray, RayPath *path, double *sLimits);
/* Find the intersections of a ray with an edge (z,r); return true if
* and only if the ray cuts the edge segment from left to right.
*/
extern int
ExitEdge(Ray *ray, double z[2], double r[2], /* inputs */
int *after, /* int *notafter, */ /* updates */
RayEdgeInfo *info); /* output */
/* Compute the path length from the current point on the ray to the
* exit (lefttoright) intersection computed by ExitEdge.
*/
extern double RayPathLength(Ray *ray, RayEdgeInfo *info);
/* Compute the path length from the exit point to the entry point (>0
* if entry AFTER exit) for the intersections computed by ExitEdge.
*/
extern double RayPathDifference(RayEdgeInfo *info);
/* Given a ray with current point entering a particular zone and side of
* the mesh, find the corresponding exit. The path length ds is >=0
* unless the ray segment lies in the negative area part of a bowtied
* zone, in which case ds<0. Return the exit side index (0<=side<=3).
* Update the current point on the ray to the exit point.
*/
extern int
ExitZone(Mesh *mesh, long zone, /* input mesh and zone number */
int side, /* input side number on entry */
Ray *ray, /* updated to exit point */
RayEdgeInfo *info[4], /* info[3] is entry side on entry,
* updated to exit side on exit */
double *ds, /* output ray path length */
double *f); /* output edge fraction */
/* Although algebraically exact formulas are used for the ExitEdge
* calculation, roundoff errors may make the redundant information
* in the Ray structure (y,z,x, and r) inconsistent. PolishExit
* attempts to restore precise selfconsistency of the ray, without
* moving the point off of the current edge (remembered in info).
*/
extern void PolishExit(Ray *ray, RayEdgeInfo *info, double *ds, double *f);
/* Through roundoff errors, it is possible for the exit point found
* in ExitZone (or FindEntryPoints) to lie slightly beyond the
* endpoints of the exit edge. AdjustRayXY moves the ray back to
* the nearby endpoint when this occurs.
*/
extern void AdjustRayXY(Ray *ray, double* z, double *r);
/* Last ditch effort by ExitZone to find the exit edge */
extern int FindLostRay(Ray *ray, RayEdgeInfo *info[4],
double z[4], double r[4], double dsx[4]);
/* Rearrange the linked list of entry points into increasing order */
extern EntryPoint *EntrySort(EntryPoint *entry);
/* Make ray path arrays longer */
extern void ExtendRayPath(RayPath *path, long pathinc);
/* Free all pointers in a RayPath (but NOT the input path pointer). */
extern void EraseRayPath(RayPath *path);
/* Adjust point on ray (x,y,z) and (z,r) to be point such that normal
* plane to ray passes through origin.
*/
void NormalizeRay(Ray *ray);
/*  */
/********** DEFAULT VALUES SET IN track.c ********/
/* Global variables control the root polishing routine PolishExit */
extern int polishRoot;
extern double polishTol1;
extern double polishTol2;
/* Global variable controlling roundoff tolerance in FindLostRay */
/* Setting this too large could result in an infinite loop, but in
* any such situation, FindLostRay will be unable to find the ray
* no matter what the setting of findRayTol */
extern double findRayTol;
/*  */
#endif
