src/snark/posit.cpp

   1 /*
   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) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  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.
  17  */
  18
  19 #include <cstdlib>
  20 #include <cstdio>
  21 #include <cmath>
  22
  23 #include "blkdta.h"
  24 #include "geom.h"
  25 #include "uiod.h"
  26
  27 #include "anglst.h"
  28 #include "posit.h"
  29
  30 void posit(INTEGER np, INTEGER nr, REAL* ax, REAL* ay, REAL* mx, REAL* my)
  31 {
  32         REAL costh;
  33         REAL sinth;
  34
  35         REAL aa, rwidth;
  36
  37         INTEGER ind;
  38
  39 // CHECK VALIDITY OF PARAMETERS NP AND NR
  40
  41         if ((np < 0) || (np >= GeoPar.prjnum))
  42         {
  43                 fprintf(output, "\n **** np = %d is out of range (0..%d)", np,
  44                                 GeoPar.prjnum - 1);
  45                 fprintf(output, "\n **** program aborted\n");
  46                 exit(777);
  47         }
  48
  49         if ((nr < 0) || (nr >= GeoPar.nrays))
  50         {
  51                 fprintf(output, "\n **** nr = %d is out of range (0..%d)", nr,
  52                                 GeoPar.nrays - 1);
  53                 fprintf(output, "\n **** program aborted\n");
  54                 exit(777);
  55         }
  56
  57         // INITIALIZATION
  58         sinth = Anglst.bsin[np];
  59         costh = Anglst.bcos[np];
  60
  61         ind = nr - GeoPar.midray;
  62
  63         // BREAK HERE FOR DIVERGENT OR PARALLEL GEOMETRY
  64
  65         if (!GeoPar.par)
  66         {
  67
  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.
  72
  73                 *ax = GeoPar.radius * costh;
  74                 *ay = GeoPar.radius * sinth;
  75
  76                 if (ind < 0)
  77                 {
  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];
  82                 }
  83                 else if (ind == 0)
  84                 {
  85                         *mx = costh;
  86                         *my = sinth;
  87                 }
  88                 else if (ind > 0)
  89                 {
  90                         *mx = costh * Anglst.cphi[ind - 1] + sinth * Anglst.sphi[ind - 1];
  91                         *my = sinth * Anglst.cphi[ind - 1] - costh * Anglst.sphi[ind - 1];
  92                 }
  93         }
  94         else
  95         {
  96
  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
 100                 // PROJECTION ANGLE.
 101
 102                 if (GeoPar.uni)
 103                         rwidth = GeoPar.pinc;
 104                 if (GeoPar.vri)
 105                         rwidth = MAX0((REAL) fabs(sinth), (REAL) fabs(costh))
 106                                         * GeoPar.pinc;
 107                 aa = rwidth * ind;
 108                 *ax = -sinth * aa;
 109                 *ay = costh * aa;
 110                 *mx = costh;
 111                 *my = sinth;
 112         }
 113
 114
 115         if (trace > 7)
 116         {
 117                 fprintf(output, "\n          posit   ax = %25.20f ay = %25.20f", *ax,
 118                                 *ay);
 119                 fprintf(output, "\n                  mx = %25.20f my = %25.20f", *mx,
 120                                 *my);
 121         }
 122         return;
 123 }