;+
;Name:
;	XRT_RT
;Purpose:
;	Calculate the spot diagram for XRT.
;Calling Sequace:
;	.run xrt_rt
;		then input the defocus value (mm)
;History:
;	1999-06-10	R.Kano	created
;-


;+
function graz_refl,a0,w0,z0,xx,uu

;Name:
;	graz_refl
:Purpose:
;	Calculate the ray reflected by a grazing mirror.
;Inputs:
;	a0:	Polynominal Coefficient of Mirror Surface
; 			r^2 = a0(0) + a0(1) * z + a0(2) * z^2 
;			            + a0(3) * z^3 + a0(4) * z^4 
;	w0:	Half Width of Mirror
;	z0:	Position of Mirror
;
;Inputs & Outputs:
;	xx:	Position Vector of Ray
;	uu:	Direction Vector of Ray
;History:
;	1999-06-10	R.Kano	created
;-
err=0

b0=a0(1:*)*[1,2,3,4]/2.

xx(2)=xx(2)-z0
xx=xx-uu*xx(2)/uu(2)

coeff=a0-[xx(0)^2+xx(1)^2,2*(uu(0)*xx(0)+uu(1)*xx(1))/uu(2),(uu(0)^2+uu(1)^2)/uu(2)^2,0,0]

z=fz_roots(coeff,/double)

s = (where((double(z) ge -w0)*(double(z) le w0)*(imaginary(z) eq 0)))(0)
if s ne -1 then begin
    z =double((z(s))(0))
    t =z/uu(2)
    xx=xx+uu*t
    nn=[xx(0:1),-poly(z,b0)] & nn=nn/sqrt(total(nn*nn))
    if total(uu*nn) gt 0 $
      then uu=uu-2*total(uu*nn)*nn $
      else err=2
endif else err=1

xx(2)=xx(2)+z0

return,err
end

;
;===== Main =====
;

; General Constants
;
am = !dpi/180 /60


; Modified Paraboloid
;
rp=170.3297306443538d0          ; mm	Radius
wp=20d0                         ; mm	Half Width
zp=-21d0                        ; mm	Position
a=[	 1.0, $                 ;	Polynominal Coefficient
	-3.6852393065221066d-3, $;	  r^2 = rp^2 * ( a0 + a1 * z'
	 9.4289065016432193d-7, $;	        + a2 * z'^2 + a3 * z'^3
	 4.2147163328109410d-7, $;	        + a4 * z'^4 )
	 1.8871992371158436d-7   ];	     z'=z/wp
ap=rp*rp*a/wp^[0,1,2,3,4]

; Modified Hyperboloid
;
rh=169.0087795268988            ; mm	Radius
wh=20d0                         ; mm	Half Width
zh=21d0                         ; mm	Position
a=[	 1.0, $                 ;	Polynominal Coefficient
	-1.1177145610964693d-2, $;	  r^2 = rp^2 * ( a0 + a1 * z'
	 2.6947219350535326d-5, $;	        + a2 * z'^2 + a3 * z'^3
	 4.2835959627179453d-7, $;	        + a4 * z'^4 )
	-1.9122005536445274d-7   ];	     z'=z/wp
ah=rh*rh*a/wh^[0,1,2,3,4]

; CCD
;
p =0.013d0                      ; mm
ans=0.0 & read,'Input the value of defocus (mm) = ',ans
zc=2700d0+ans                   ; mm


;z=findgen(4001)/100-20
;plot ,z-21,sqrt(poly(z,ap)) $
;  ,yra=[-175,175],xra=[-50,3000]
;;  ,yra=[160,175],xra=[-50,50]
;oplot,z+21,sqrt(poly(z,ah))



!p.multi=[0,3,3]
!p.charsize=2
for i=0,16,2 do begin

    spt=[0,0]
    err=1

    dr=0.05d0
    for r=170.0d0,170.7d0,dr do begin

        for phi=0d0,2*!dpi,20*2*!dpi/(r/dr) do begin

            xx=[r*cos(phi),r*sin(phi),-50]
            uu=[0,sin(i*am),cos(i*am)]

            xx0=xx	; Initial Ray
            uu0=uu
            ep =graz_refl(ap,wp,zp,xx,uu)
            xx1=xx	; Ray Reflected by Paraboloid
            uu1=uu
            eh =graz_refl(ah,wh,zh,xx,uu)
            xx2=xx	; Ray Reflected by Hyperboloid
            uu2=uu

            xx3=xx2-uu2*(xx2(2)-zc)/uu2(2)	; Coordinates at CCD
            spt=[[spt],[xx3(0:1)]]
            err=[[err],[ep+10*eh]]
            
;           oplot,[xx0(2),xx1(2)],[xx0(0),xx1(0)]
;           if ep eq 0 then begin
;               oplot,[xx1(2),xx2(2)],[xx1(0),xx2(0)]
;               if eh eq 0 $
;                   then oplot,xx2(2)+uu2(2)*[0,1e4],xx2(0)+uu2(0)*[0,1e4] $
;                   else plots,xx2(2),xx2(0),psym=7
;           endif else plots,xx1(2),xx1(0),psym=7
        
        endfor
    endfor
    
    ss=where(err eq 0)
    ns=n_elements(ss)
    ms=total(spt(*,ss),2)/ns
    rs=sqrt(total(spt(*,ss)^2)/ns-total(ms^2))
    text1=string(ms,form='(F6.2," mm/",F6.2," mm")')
    text2=string(2*rs*1000,form='(F5.1," um =")') $
      +string(2*rs/0.013,form='(F5.2," arcsec")')
    print,text1,' : ',text2

    plot,p*[-1,-1,1,1,-1]/2,p*[-1,1,1,-1,-1]/2 $
      ,xra=3*p*[-1,1],yra=3*p*[-1,1] $
      ,xmargin=[10,7] $
      ,title='offpoint='+strtrim(i,2)+'arcmin',xtit=text1,subtit=text2
    oplot,spt(0,ss)-ms(0),spt(1,ss)-ms(1),psym=3

endfor

end