; pro stardust_geometry, sc_sun_pos,sc_target_pos,ra,dec,twist,$ ; sunca,nepca,ncpca ; ; program to compute the projected sunward direction and the ; projected North pole directions for stardust images ; ; Input Parameters: ; ; sc_sun_pos[3] real from the image header ; sc_target_pos[3] real from the image header ; ra real from the image header ; dec real from the image header ; twist real from the image header ; ; Output Parameters: ; sunca real sunward clock angle ; nepca real North Ecliptic pole clock angle ; ncpca real North Celestial pole clock angle ; ;==================================================================== function axisrot,vect,iaxis,angle ; function to rotate a set of coordinates about one of the axes ; Note: configured so that rotation is by the right-hand rule ; definition around the given axis ; ; Input: ; vect(3) real three coordinates to be transformed ; iaxis integer axis to be rotated about (0,1,2) ; angle real angle to rotate (radians) ; ; Output: ; vect2(3) real three coordinates in the new frame on_error,2 vect2 = fltarr(3) ca = cos(angle) sa = -1. * sin(angle) case iaxis of 0: begin vect2(0) = vect(0) vect2(1) = vect(1)*ca - vect(2)*sa vect2(2) = vect(1)*sa + vect(2)*ca end 1: begin vect2(0) = vect(2)*sa + vect(0)*ca vect2(1) = vect(1) vect2(2) = vect(2)*ca - vect(0)*sa end 2: begin vect2(0) = vect(0)*ca - vect(1)*sa vect2(1) = vect(0)*sa + vect(1)*ca vect2(2) = vect(2) end endcase return,vect2 end ;==================================================================== pro stardust_geometry, sc_sun_pos,sc_target_pos,ra,dec,twist,$ sunca,nepca,ncpca ; compute the celestial north clock angle direction r2d = 180./!pi cel_n_clkang = 180. - twist if cel_n_clkang lt 0. then cel_n_clkang = cel_n_clkang+360. ; compute the vector from comet to sun target_sun_pos = sc_sun_pos - sc_target_pos ; set up rotation angles rotangr = !pi/2. + ra/r2d rotangd = !pi/2. - dec/r2d rotangt = (180.-twist)/r2d ; rotate sun vector to plane of the sky reference frame vec1 = target_sun_pos vec2 = axisrot(vec1,2,rotangr) vec1 = axisrot(vec2,0,rotangd) vec2 = axisrot(vec1,2,-rotangt) ; compute phase angle as a check r2 = sqrt(vec2[0]^2 +vec2[1]^2 +vec2[2]^2 ) phase = !pi - acos(vec2[2]/r2) ; compute the sunward clock angle proj_sun_clkang = atan(-vec2[0],vec2[1]) if proj_sun_clkang lt 0. then proj_sun_clkang = proj_sun_clkang+2*!pi ; compute the ecliptic North vector obliquity = 23.46 ; earth's obliquity at about the time of encounter ncp = [0.,0.,1.] nep = axisrot(ncp,0,-obliquity/r2d) ; rotate NEP to the plane of the sky coordinates vec1 = nep vec2 = axisrot(vec1,2,rotangr) vec1 = axisrot(vec2,0,rotangd) vec2 = axisrot(vec1,2,-rotangt) ; compute the NEP clock angle proj_nep_clkang = atan(-vec2[0],vec2[1]) if proj_nep_clkang lt 0. then proj_nep_clkang = proj_nep_clkang+2*!pi f1='(a,f7.2)' print,' ' print,'Projected Solar Vector Clock Angle: ',proj_sun_clkang*r2d,format=f1 print,'Solar Phase Angle: ',phase*r2d,format=f1 print,'North Ecliptic Pole Clock Angle: ',proj_nep_clkang*r2d,format=f1 print,'North Celestial Pole Clock Angle: ',cel_n_clkang,format=f1 print,' ' sunca = proj_sun_clkang*r2d nepca = proj_nep_clkang*r2d ncpca = cel_n_clkang return end