pro for_global_local_rotate,chromoheight,theta_prime,phi_prime,r3D,thmod,phmod,theta_global,phi_global

; 
; Name: FOR_GLOBAL_LOCAL_ROTATE
;
;  rotates from local vertical (x',y',z) to global 3D (x,y,z) coordinate system
;
; INPUTS:
;
;  chromoheight, theta_prime, phi_prime are the local spherical coordinates of cone footpoints
;           (note called theta_prime_2D, phi_prime_2D in for_uv_stokes.pro)
;
;  r3D, thmod, phmod are the global spherical coordinates of the scattering point P 
;       **IN MODEL FRAME**
;	 *(rotated by B angle and central meridian and potentially thetao)*
;
; OUTPUTS
;    theta_global, phi_global are the global **MODEL FRAME, not observer** spherical coordinate 
;           positions of cone footpoints
;  (note called theta3d_chromo, phi3d_chromo in for_uv_stokes.pro)
;    (radial position is still chromoheight)
;
; Called by FOR_UV_STOKES
;
; Based on notes from Roberto Casini
;  *BUT using FORWARD convention where primed coordinates are in local vertical frame
;    and made to be consistent with local-global relationship as defined in FOR_UVANGLE.PRO
;
        st=sin(theta_prime)
        ct=cos(theta_prime)
        sp=sin(phi_prime)
        cp=cos(phi_prime)

; this is in local (primed) coordinate system 
;  representing footpoints of the cone starting from scattering point P
;  _C stands for cone  
	xp_C = chromoheight*st*cp
        yp_C = chromoheight*st*sp 
        zp_C = chromoheight*ct 

; this is in global coordinate system
;  representing Cartesian coordinates of scattering point P
;  _P stands for scattering point P
	x_P = r3D*sin(thmod)*cos(phmod)
	y_P = r3D*sin(thmod)*sin(phmod)
	z_P = r3D*cos(thmod)

; these are equivalent to LHS of final expression of eq. 16 in Roberto's notes

; following local-global rotation defined in FOR_UVANGLE.PRO
;
	alpha = -1.d0*atan(x_P,sqrt(y_P^2+z_P^2))
 	beta = atan(y_P,z_P)

; Now rotating the cone footpoints backward from the rotation done in FOR_UVANGLE.PRO
; i.e., from local to global so taking the inverse of R_SS'=inv(matrix(alpha)*matrix(beta)) from eq. 1
; NOTE not yet translated 
        x_C_nt = cos(double(alpha))*xp_C - sin(double(alpha))*zp_C

        y_C_nt = sin(double(alpha))*sin(double(beta))*xp_C + cos(double(beta))*yp_C + cos(double(alpha))*sin(double(beta))*zp_C

        z_C_nt = sin(double(alpha))*cos(double(beta))*xp_C - sin(double(beta))*yp_C + cos(double(alpha))*cos(double(beta))*zp_C

; these are equivalent to square bracket in RHS of eq. 16
; 
        xi = r3D/chromoheight

	eq_17=theta_prime*0.d0

        testtp=where(theta_prime ne 0, ntp)
	if ntp ne 0 then begin
          cot2_th = 1./(tan(theta_prime[testtp]))^2
	  cos_delS = (xi[testtp] + sqrt((cot2_th*(1. + cot2_th - xi[testtp]^2))))/(1. + cot2_th)
	  eq_17[testtp] = (xi[testtp] - cos_delS)/cos(theta_prime[testtp])
        endif 

        testtpz=where(theta_prime eq 0, ntpz)
        if ntpz ne 0 then begin
; special case points lying along local vertical
;   theta_prime=0, cos_delS=1
          eq_17[testtpz] = xi[testtpz] - 1.d0
        endif

; Now solving eq.16 for translated cone footpoints in Cartesian global coordinates
	x_C = x_P - eq_17 * x_C_nt

	y_C = y_P - eq_17 * y_C_nt
	
	z_C = z_P - eq_17 * z_C_nt

;  this next should be equivalent
;
     ;   rat=z_C/chromoheight
     ;   testpos = where(rat gt 1.)
     ;   testneg = where(rat lt -1.)
     ;   if min(testpos) ne -1 then rat[testpos]=1.
     ;   if min(testneg) ne -1 then rat[testneg]=-1.

     ;   theta_global = acos(rat)
        theta_global = atan(sqrt(x_C^2+y_C^2),z_C)

        phi_global = atan(y_C,x_C)

; this next for debugging
        test = where(finite(theta_global) eq 0.,ntt)
        if ntt ne 0. then begin
         print,'warning! NaN in theta_global' 
         stop
        endif
        test = where(finite(phi_global) eq 0.,ntp)
        if ntp ne 0. then begin
          print,'warning! NaN in phi_global' 
          stop
        endif
end