2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/posit.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 POSIT IS USED TO DETERMINE THE EQUATION OF A RAY. ITS INPUT
11 CONSISTS OF THE NUMBER OF THE RAY AND THE SINE AND COSINE OF
12 THE PROJECTION ANGLE. IT RETURNS THE COORDINATES OF A POINT
13 ON THE RAY, (AX,AY), AND THE SLOPE OF THE RAY, (MX,MY). THE
14 ROUTINE ACCESSES GEOMETRY INFORMATION CONTAINED IN /GEOM/ AND
15 ASSUMES THAT GENPHI HAS BEEN CALLED TO COMPUTE THE SINE AND COSINE
16 TABLES FOR THE DIVERGENT GEOMETRY CASE.
30 void posit(INTEGER np, INTEGER nr, REAL* ax, REAL* ay, REAL* mx, REAL* my)
39 // CHECK VALIDITY OF PARAMETERS NP AND NR
41 if ((np < 0) || (np >= GeoPar.prjnum))
43 fprintf(output, "\n **** np = %d is out of range (0..%d)", np,
45 fprintf(output, "\n **** program aborted\n");
49 if ((nr < 0) || (nr >= GeoPar.nrays))
51 fprintf(output, "\n **** nr = %d is out of range (0..%d)", nr,
53 fprintf(output, "\n **** program aborted\n");
58 sinth = Anglst.bsin[np];
59 costh = Anglst.bcos[np];
61 ind = nr - GeoPar.midray;
63 // BREAK HERE FOR DIVERGENT OR PARALLEL GEOMETRY
68 // DIVERGENT CALCULATIONS FOLLOW
69 // RETURN THE POSITION OF THE SOURCE AS THE POINT ON THE RAY.
70 // THE ANGLE OF THE RAY IS THE DIFFERENCE OF THE PROJECTION ANGLE
71 // AND THE ANGLE THE RAY MAKES WITH THE CENTER RAY.
73 *ax = GeoPar.radius * costh;
74 *ay = GeoPar.radius * sinth;
78 *mx = costh * Anglst.cphi[(-ind) - 1]
79 - sinth * Anglst.sphi[(-ind) - 1];
80 *my = sinth * Anglst.cphi[(-ind) - 1]
81 + costh * Anglst.sphi[(-ind) - 1];
90 *mx = costh * Anglst.cphi[ind - 1] + sinth * Anglst.sphi[ind - 1];
91 *my = sinth * Anglst.cphi[ind - 1] - costh * Anglst.sphi[ind - 1];
97 // PARALLEL CALCULATIONS FOLLOW
98 // COMPUTE THE INTERSECTION OF THE RAY WITH A LINE PERPENDICULAR TO
99 // IT AND PASSING THROUGH THE ORIGIN. THE ANGLE IS JUST THE
103 rwidth = GeoPar.pinc;
105 rwidth = MAX0((REAL) fabs(sinth), (REAL) fabs(costh))
117 fprintf(output, "\n posit ax = %25.20f ay = %25.20f", *ax,
119 fprintf(output, "\n mx = %25.20f my = %25.20f", *mx,