# NAIF user

Linda Wills (linda@cc.gatech.edu)
28 Oct 1994 13:37:40 -0400

Hi! Here's a reply from the NAIF client in San Diego. -- Linda

~From: mc@msss.com (Mike Caplinger)
To: Linda Wills <linda@ee.gatech.edu>
~Subject: Re: Donating FORTRAN code to science
~Date: Wed, 26 Oct 94 17:07:59 PDT

Linda,

We have a fairly significant collection of NAIF-using code, but it's
in C, not Fortran -- we developed a bunch of interface macros that
allow us to use NAIF conveniently from C using the Unix Fortran
calling conventions.

Here's a piece of code in Fortran that represents the average
complexity of what we do. No warranty, etc. This is old code; I
think we do things more cleanly now. Let me know if there's anything
you can do with the C code.

subroutine spice

integer bodyid

double precision lat1, lon1
double precision xsc(6), tipm(3, 3), subsc(3), axes(3)
double precision normal(3)
double precision cang(3), pang(3)
double precision a, b, c, rho, alt, lontmp, lattmp
double precision camvec(3)
double precision cpoint(3)
double precision halfpi, dpr, rpd
double precision range, ra, dec, w
double precision cmat(3,3)
double precision big(3,3)
double precision xscin(3), cangin(3), pangin(3)
logical found
double precision try(3)
double precision fl
integer x, y, cx, cy

return

entry spinit(xscin, cangin, pangin)

fl = 40282.0
cx = 1204/2
cy = 1056/2

bodyid = 4

axes(1) = 3393.4
axes(2) = 3393.4
axes(3) = 3375.8

a = axes(1)
b = axes(2)
c = axes(3)

do i = 1,3
xsc(i) = xscin(i)
cang(i) = cangin(i)
pang(i) = pangin(i)
enddo

c convert camera angles to rectangular coords and normalize
call sphrec(1d3, halfpi() - cang(1), cang(2), camvec)
call unorm(camvec, camvec, w)
c print *, 'raw cam ', camvec

call rotate(pang(2)+halfpi(), 3, tipm)
call rotmat(tipm, halfpi()-pang(1), 1, tipm)
call rotmat(tipm, pang(3), 3, tipm)

call mxv(tipm, xsc, xsc)
call mxv(tipm, camvec, camvec)

c build matrix to translate from S/C to camera coord system...
c right ascension first
call rotate(halfpi()+cang(2), 3, cmat)
c then declination
call rotmat(cmat, halfpi()-cang(1), 1, cmat)
c finally twist
call rotmat(cmat, cang(3), 3, cmat)

c big takes a vector from camera coords to body-fixed coords
c in one big jump, and its transpose takes it back.
call mxmt(tipm, cmat, big)

return

entry spcenter(lon1, lat1)

c calculate image center coordinates
call surfpt(xsc, camvec, a, b, c, cpoint, found)
if(found) then
c print *, 'found intersection'
call reclat(cpoint, rho, lon1, lat1)
call surfnm(a, b, c, cpoint, normal)
call reclat(normal, rho, lontmp, lat1)
lat1 = lat1 * dpr()
lon1 = lon1 * dpr()
lon1 = -lon1
endif

return

entry splonlat(x, y, lon1, lat1)

c big now takes a vector from camera coords to body-fixed coords in one big jump...
c call mxv(big, pyr, pyr)
c ...and its transpose takes it back.
c call mtxv(big, pyr, pyr)

camvec(3) = fl
camvec(1) = x-cx
camvec(2) = y-cy
call vhat(camvec, camvec)

call mxv(big, camvec, camvec)
call surfpt(xsc, camvec, a, b, c, cpoint, found)
if(found) then
c print *, 'found intersection'
call reclat(cpoint, rho, lon1, lat1)
call surfnm(a, b, c, cpoint, normal)
call reclat(normal, rho, lontmp, lat1)
lat1 = lat1 * dpr()
lon1 = lon1 * dpr()
lon1 = -lon1
if(lon1 .lt. 0.0) lon1 = lon1 + 360.0
endif

return

entry spxy(lon1, lat1, x, y)

lon1 = -lon1
lat1 = lat1 * rpd()
lon1 = lon1 * rpd()

call georec(lon1, lat1, 0.0, axes(1),
+(axes(1)-axes(3))/axes(1), cpoint)

call vsub(cpoint, xsc, try)

c need to test this to see if it's on the right side of the planet!!!
call mtxv(big, try, try)
alpha = try(3)/fl
try(1) = try(1)/alpha
try(2) = try(2)/alpha

x = try(1) + cx + 0.5
y = try(2) + cy + 0.5
end

```--

-Linda Wills (linda@cc.gatech.edu)
```