~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 degrees to radians
call convrt(cang(1), 'DEGREES', 'RADIANS', cang(1))
call convrt(cang(2), 'DEGREES', 'RADIANS', cang(2))
call convrt(cang(3), 'DEGREES', 'RADIANS', cang(3))
call convrt(pang(1), 'DEGREES', 'RADIANS', pang(1))
call convrt(pang(2), 'DEGREES', 'RADIANS', pang(2))
call convrt(pang(3), 'DEGREES', 'RADIANS', pang(3))
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)