function sphtrig_translb,theta,phi,tilt=tilt,twist=twist,eps=eps,\$ verbose=verbose, _extra=e ;+ ;function sphtrig_translb ; transform the given {theta,phi} coordinates on a sphere that ; has been tilted and twisted (in that order) to {theta',phi'} ; on an untwisted and untilted pov reference frame, and ; returns them in a structure. ; ;syntax ; tpstr=sphtrig_translb(theta,phi,tilt=tilt,twist=twist,\$ ; eps=eps,verbose=vv) ; ;parameters ; theta [INPUT; required] angular coordinate from Z-axis [deg] ; phi [INPUT; required] azimuthal coordinate around Z-axis [deg] ; * THETA,PHI may be vectors or scalars ; * if vectors, assumed to be paired sets, with ; values filled in with 0's if one falls short ; ;keywords ; tilt [INPUT; default=0] inclination angle of sphere ; relative to pov coordinate frame [deg] ; twist [INPUT; default=0] rotation angle by which the ; vector going through the pole is rotated around ; the Z-axis of the pov reference frame [deg] ; * note that TWIST is given in the negative sense -- ; the angle by which the frame must be rotated to ; return to 0. the program takes the -ve of the ; supplied TWIST. ; eps [INPUT; default=1e-6] a small number ; * There are some singularities when THETA and PHI are zero. ; This is added on to make sure at least something ; close to the limiting value is returned in that case ; verbose [INPUT] controls chatter ; _extra [JUNK] here only to prevent crashing the program ; ;example ; .run sphtrig_translb ; makes lines of constant latitude, constant longitude, and ; draws the figure used to derive the transformation used here. ; (the ultimate in self-documenting!) ; ;subroutines ; example run requires PEASECOLR ; ;history ; vinay kashyap (2013feb) ; added keyword EPS (2013mar) ; corrected sign error that was making PHI go clockwise instead of ; counter-clockwise (2015sep) ;- ; usage ok='ok' & npar=n_params() & nt=n_elements(theta) & np=n_elements(phi) if npar lt 2 then ok='Insufficient parameters' else \$ if nt eq 0 then ok='THETA not defined' else \$ if np eq 0 then ok='PHI not defined' if ok ne 'ok' then begin print,'Usage: tpstr=sphtrig_translb(theta,phi,tilt=tilt,twist=twist,eps=eps,verbose=vv)' print,' compute theta,phi transformed from tilted and twisted sphere' print,' to line up with point-of-view coordinate frame' if npar gt 0 then message,ok,/informational return,-1L endif ; parse input nn=nt > np zth=dblarr(nn) & zth[0L:np-1L]=theta[*] zph=dblarr(nn) & zph[0L:nt-1L]=phi[*] ; parse keywords vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L if not keyword_set(eps) then eps=1d-6 tht0=0.D + eps & if keyword_set(tilt) then tht0=tilt[0] while tht0 lt 0 do tht0=tht0+180.D if tht0 gt 180. then tht0=(tht0 mod 180.D) phi0=0.D + eps & if keyword_set(twist) then phi0=twist[0] while phi0 lt 0 do phi0=phi0+360.D if phi0 gt 360. then phi0=(phi0 mod 360.D) phi0=-phi0 ; convert to radians zth=zth*!dpi/180.D & zph=zph*!dpi/180.D & tht0=tht0*!dpi/180.D & phi0=phi0*!dpi/180.D ; points, arcs, and angles refer to the documentation figure, which ; can be viewed by running the example script to drive this program. ; first, get angle CSB ;arcBC=2.*!DPI-zph ;this changes sign of phi arcBC=zph arcSC=!DPI/2. arcSB=!DPI/2. ; cos(arcBC) = cos(arcSC)*cos(arcSB) + sin(arcSC)*sin(arcSB)*cos(angCSB) ; ==> cos(arcBC) = 0*0 + 1*1*cos(angCSB) ; ==> angCSB = arcBC angCSB=arcBC ; now, get arc ZX arcSX=zth arcZS=tht0 angZSX=!pi-angCSB ; cos(arcZX) = cos(arcZS)*cos(arcSX) + sin(arcZS)*sin(arcSX)*cos(angZSX) arcZX = acos(cos(arcZS)*cos(arcSX) + sin(arcZS)*sin(arcSX)*cos(angZSX)) sgnZSO = intarr(nn)+1 & o0=where(angZSX lt 0,mo0) & if mo0 gt 0 then sgnZSO[o0]=-1 ; get angle XZS ; cos(arcSX) = cos(arcZS)*cos(arcZX) + sin(arcZS)*sin(arcZX)*cos(angXZS) ; cos(angXZS) = (cos(arcSX)-cos(arcZS)*cos(arcZX))/(sin(arcZS)*sin(arcZX)) angXZS = sgnZSO*acos((cos(arcSX)-cos(arcZS)*cos(arcZX))/(sin(arcZS)*sin(arcZX))) ; angle UZV arcZU=!DPI/2. arcZV=!DPI/2. arcUV=phi0 ; cos(arcUV) = cos(arcZU)*cos(arcZV) + sin(arcZU)*sin(arcZV)*cos(angUZV) ; ==> cos(arcUV) = 0*0 + 1*1*cos(angUZV) ; ==> angUZV = arcUV angUZV=arcUV ; angle WZU angWZU = angXZS - angUZV ; arc UW arcZW=!DPI/2. ; cos(arcUW) = cos(arcZW)*cos(arcZU) + sin(arcZW)*sin(arcZU)*cos(angWZU) ; ==> cos(arcUW) = 0*0 + 1*1*cos(angWZU) arcUW=angWZU ; output ; convert to [deg] tpstr=create_struct('THETA',arcZX*180./!pi,'PHI',arcUW*180./!pi) if vv gt 1000 then stop,'HALTing; type .CON to continue' return,tpstr end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; example ; set variables TILT, TWIST, STHETA, SPHI ; ; made slightly more verbose (VK; 2015sep) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; peasecolr & loadct,3 & peasecolr ; calling sequence print,'' & print,'CALLING SEQUENCE:' & print,'' jnk=sphtrig_translb() print,'' ; initialize if n_elements(tilt) ne 1 then tilt=35. else tilt=tilt[0] if n_elements(twist) ne 1 then twist=50. else twist=twist[0] if not keyword_set(verbose) then verbose=1 print,'SET THE INCLINATION OF THE STAR USING VARIABLES --' help,tilt,twist if tilt lt 0 then begin message,'what does negative TILT mean anyway?',/informational tilt=abs(tilt) help,tilt,twist endif if tilt gt 180 then begin message,'nothing can TILT beyond a vertical dive',/informational tilt=(tilt mod 180.) help,tilt,twist endif wait,1 print,'FIRST DISPLAY A LATITUDE GRID' if !d.name eq 'X' then window,0,xsize=800,ysize=800 ; set up observer's frame surface,fltarr(2,2),/nodata,xrange=1.2*[-1,1],yrange=1.2*[-1,1],zrange=1.2*[-1,1],xstyle=5,ystyle=5,zstyle=5,ax=20,az=-100,/save plots,[0,1],[0,0],[0,0],col=3,/t3d,line=1 plots,[0,0],[0,1],[0,0],col=3,/t3d,line=1 plots,[0,0],[0,0],[0,1],col=3,/t3d,line=1 ; now set up the source sphere frame plots,[0,sin(tilt*!pi/180.)*cos(twist*!pi/180.)],[0,sin(tilt*!pi/180.)*sin(twist*!pi/180.)],[0,cos(tilt*!pi/180.)],col=1,/t3d,thick=2 ; draw latitudes ang2=findgen(361) for ilat=0,180,10 do begin ang1=0*ang2+ilat tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,/t3d,psym=3,col=1 ox=where(xx gt 0,mox) & if mox gt 0 then plots,xx[ox],yy[ox],zz[ox],/t3d,psym=4,col=1 if ilat eq 90 then plots,xx,yy,zz,/t3d,col=1,thick=2 endfor wait,1 print,'NEXT DISPLAY A LONGITUDE GRID' if !d.name eq 'X' then window,2,xsize=800,ysize=800 ; set up observer's frame surface,fltarr(2,2),/nodata,xrange=1.2*[-1,1],yrange=1.2*[-1,1],zrange=1.2*[-1,1],xstyle=5,ystyle=5,zstyle=5,ax=20,az=-100,/save plots,[0,1],[0,0],[0,0],col=3,/t3d,line=1 plots,[0,0],[0,1],[0,0],col=3,/t3d,line=1 plots,[0,0],[0,0],[0,1],col=3,/t3d,line=1 ; now set up the source sphere frame plots,[0,sin(tilt*!pi/180.)*cos(twist*!pi/180.)],[0,sin(tilt*!pi/180.)*sin(twist*!pi/180.)],[0,cos(tilt*!pi/180.)],col=1,/t3d,thick=2 ; draw longitudes ang1=findgen(181) for ilon=0,340,20 do begin ang2=0*ang1+ilon tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,/t3d,psym=3,col=1 ox=where(xx gt 0,mox) & if mox gt 0 then plots,xx[ox],yy[ox],zz[ox],/t3d,psym=4,col=1 if ilon eq 0 then plots,xx,yy,zz,/t3d,col=1,thick=2 endfor ; pick a point, any point, on the source sphere if n_elements(stheta) ne 1 then stheta=44. if n_elements(sphi) ne 1 then sphi=111. while stheta lt 0 do stheta=stheta+180. if stheta gt 180. then stheta=(stheta mod 180.) while sphi lt 0 do sphi=sphi+360. if sphi gt 360. then sphi=(sphi mod 360.) wait,1 print,'SELECT AN ARBITRARY POINT ON THE STAR USING VARIABLES --' help,stheta,sphi print,'SHOW WHERE THE POINT LIES, ALONG WITH EXPLICATORY SPHERICAL TRIANGLES --' if !d.name eq 'X' then window,1,xsize=800,ysize=800 ; set up the observer's frame surface,fltarr(2,2),/nodata,xrange=1.2*[-1,1],yrange=1.2*[-1,1],zrange=1.2*[-1,1],xstyle=5,ystyle=5,zstyle=5,ax=20,az=-100,/save plots,[0,1],[0,0],[0,0],col=3,/t3d,line=1 plots,[0,0],[0,1],[0,0],col=3,/t3d,line=1 plots,[0,0],[0,0],[0,1],col=3,/t3d,line=1 xyouts,0.05,0.05,z=0.05,text_axes=2,charthick=2,charsize=2,col=3,/t3d,'O' xyouts,1.05,0,z=-0.05,text_axes=2,charthick=2,charsize=2,col=3,/t3d,'U' xyouts,0.05,0,z=1.05,text_axes=2,charthick=2,charsize=2,col=3,/t3d,'Z' ; draw some orienting circles ang=findgen(360)*!pi/180. xx=cos(ang) & yy=sin(ang) & zz=0*ang & plots,xx,yy,zz,col=3,line=1,/t3d xx=0*ang & yy=sin(ang) & zz=cos(ang) & plots,xx,yy,zz,col=3,line=1,/t3d ; mark the tilt ang=[findgen(fix(tilt)+1),tilt]*!pi/180. & xx=sin(ang)*cos(twist*!pi/180.) & yy=sin(ang)*sin(twist*!pi/180.) & zz=cos(ang) plots,xx,yy,zz,col=3,/t3d nang=n_elements(ang) & xyouts,xx[nang-1],yy[nang-1],z=zz[nang-1]+0.05,text_axes=2,charthick=2,charsize=2,col=1,/t3d,'S' ; now set up the source sphere frame plots,[0,sin(tilt*!pi/180.)*cos(twist*!pi/180.)],[0,sin(tilt*!pi/180.)*sin(twist*!pi/180.)],[0,cos(tilt*!pi/180.)],col=1,/t3d ; draw some orienting circles in source frame ang2=findgen(360) & ang1=0*ang2+90 tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist,verbose=verbose) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,psym=3,col=1,/t3d ang1=findgen(91) & ang2=0*ang1 tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist,verbose=verbose) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,line=2,col=1,/t3d nang=n_elements(ang1) & xyouts,xx[nang-1]-0.05,yy[nang-1]-0.05,z=zz[nang-1]-0.1,text_axes=2,charthick=2,charsize=2,col=1,/t3d,'B' jnk=min(abs(zz),imn) & plots,[0,xx[imn]],[0,yy[imn]],[0,zz[imn]],line=1,col=1,/t3d xyouts,xx[imn]-0.05,yy[imn]-0.05,z=zz[imn]-0.1,text_axes=2,charthick=2,charsize=2,col=1,/t3d,'V' plots,[0,xx[nang-1]],[0,yy[nang-1]],[0,zz[nang-1]],line=1,col=2,/t3d ; mark the point slb=sphtrig_translb(stheta,sphi,tilt=tilt,twist=twist,verbose=verbose) szz=cos(slb.THETA*!pi/180.) & sxx=sin(slb.THETA*!pi/180.)*cos(slb.PHI*!pi/180.) & syy=sin(slb.THETA*!pi/180.)*sin(slb.PHI*!pi/180.) plots,[0,sxx],[0,syy],[0,szz],col=2,/t3d,symsize=2,thick=2 xyouts,sxx-0.05,syy-0.05,z=szz+0.05,text_axes=2,charthick=2,charsize=2,col=2,/t3d,'X' print,'Cartesian (x,y,z) coordinates of X in upright frame:',sxx[0],syy[0],szz[0] print,'Polar (theta,phi) coordinates of X in upright frame:',(slb.THETA)[0],(slb.PHI)[0] ; mark the triangle in the source frame ang1=findgen(91) & ang2=0*ang1+sphi tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist,verbose=verbose) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,line=2,col=1,/t3d nang=n_elements(ang1) & xyouts,xx[nang-1]-0.05,yy[nang-1]-0.05,z=zz[nang-1]-0.1,text_axes=2,charthick=2,charsize=2,col=1,/t3d,'C' plots,[0,xx[nang-1]],[0,yy[nang-1]],[0,zz[nang-1]],line=1,col=2,/t3d ; ang1=[findgen(fix(stheta)+1),stheta] & ang2=0*ang1+sphi tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist,verbose=verbose) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,col=1,/t3d ; ang2=[findgen(fix(sphi)+1),sphi] & ang1=0*ang2+90. tpstr=sphtrig_translb(ang1,ang2,tilt=tilt,twist=twist,verbose=verbose) zz=cos(tpstr.THETA*!pi/180.) & xx=sin(tpstr.THETA*!pi/180.)*cos(tpstr.PHI*!pi/180.) & yy=sin(tpstr.THETA*!pi/180.)*sin(tpstr.PHI*!pi/180.) plots,xx,yy,zz,line=2,col=1,/t3d ; mark the triangle in the observer's frame ang1=findgen(91) & ang2=0*ang1+(slb.PHI)[0] zz=cos(ang1*!pi/180.) & xx=sin(ang1*!pi/180.)*cos(ang2*!pi/180.) & yy=sin(ang1*!pi/180.)*sin(ang2*!pi/180.) plots,xx,yy,zz,line=2,col=4,/t3d nang=n_elements(ang1) & xyouts,xx[nang-1]-0.05,yy[nang-1]-0.05,z=zz[nang-1]-0.1,text_axes=2,charthick=2,charsize=2,col=4,/t3d,'W' plots,[0,xx[nang-1]],[0,yy[nang-1]],[0,0],line=1,col=4,/t3d ; ang1=[findgen(fix((slb.THETA)[0])+1),(slb.THETA)[0]] & ang2=0*ang1+(slb.PHI)[0] zz=cos(ang1*!pi/180.) & xx=sin(ang1*!pi/180.)*cos(ang2*!pi/180.) & yy=sin(ang1*!pi/180.)*sin(ang2*!pi/180.) plots,xx,yy,zz,col=4,/t3d ; ang1=findgen(91) & ang2=0*ang1 zz=cos(ang1*!pi/180.) & xx=sin(ang1*!pi/180.)*cos(ang2*!pi/180.) & yy=sin(ang1*!pi/180.)*sin(ang2*!pi/180.) plots,xx,yy,zz,line=2,col=4,/t3d ; some descriptive labels xyouts,0.05,0.95-0.*0.03,/norm,charthick=1.5,charsize=2,'ZS = TILT = '+string(tilt,'(f7.2)') xyouts,0.05,0.95-1.*0.03,/norm,charthick=1.5,charsize=2,'UOV = TWIST = '+string(twist,'(f7.2)') xyouts,0.05,0.95-2.*0.03,/norm,charthick=1.5,charsize=2,'SX = STHETA = '+string(stheta,'(f7.2)') xyouts,0.05,0.95-3.*0.03,/norm,charthick=1.5,charsize=2,'COB = SPHI = '+string(sphi,'(f7.2)') xyouts,0.05,0.95-5.*0.03,/norm,charthick=1.5,charsize=2,'WX = LAT = '+string(90.-(slb.theta)[0],'(f7.2)') xyouts,0.05,0.95-6.*0.03,/norm,charthick=1.5,charsize=2,'UOW = LON = '+string((slb.phi)[0],'(f7.2)') end