add_error_on_gp.pro (Kronos IDL library)

;+
; Contains the add_error_on_gp procedure
;
; :Author:
;    Baptiste Cecconi, Laurent Lamy
;
; :History:
;    2008/03/05: Created
;
;    2008/03/05: Last Edit
;-
;
;+
; add_error_on_gp is a procedure that <behavior desc here>
;
; :Params:
;    dth:    in, required, type=sometype
;       A parameter named dth
;    rr:     in, required, type=sometype
;       A parameter named rr
;    ydf:    in, required, type=sometype
;       A parameter named ydf
;    zdf:    in, required, type=sometype
;       A parameter named zdf
;    qrot:   in, required, type=sometype
;       A parameter named qrot
;    nerr:   in, required, type=sometype
;       A parameter named nerr
;    qdferr: in, required, type=sometype
;       A parameter named qdferr
;
; :Keywords:
;    axis:   in, optional, type=sometype
;       A keyword named axis
;-
pro add_error_on_gp,dth,rr,ydf,zdf,qrot,nerr,qdferr,axis=axis

ndata = max([n_elements(dth),n_elements(rr),n_elements(ydf),n_elements(zdf)])

qdferr = fltarr(4,nerr,ndata)
zdferr = fltarr(nerr,ndata)
ydferr = fltarr(nerr,ndata)
azierr = fltarr(nerr,ndata)
laterr = fltarr(nerr,ndata)

r0=rr*tan(dth)

if ~ keyword_set(axis) then axis='N'
case axis of 
  'N'  : phi0 = !pi/2.
  'NW' : phi0 = 3.*!pi/4.
  'W'  : phi0 = !pi
  'SW' : phi0 = 5.*!pi/4.
  'S'  : phi0 = 3.*!pi/2.
  'SE' : phi0 = 7.*!pi/4.
  'E'  : phi0 = 0.
  'NE' : phi0 = !pi/4.
end

for i=0,nerr-1 do begin
  zdferr(i,*) = reform(zdf + sin(float(i)/nerr*2*!pi+phi0)*r0,1,ndata)
  ydferr(i,*) = reform(ydf + cos(float(i)/nerr*2*!pi+phi0)*r0,1,ndata)

  azierr(i,*) = reform(atan(-ydferr(i,*),rr),1,ndata)      
  laterr(i,*) = reform(atan(zdferr(i,*),rr),1,ndata)

;  qdferr([1:3],i,*) = reform(make_vect_sph(1.,reform(!pi/2.-laterr(i,*),ndata),reform(azierr(i,*),ndata)),3,1,ndata)
  qdferr(3,i,*) =  reform(sin(laterr(i,*)),1,1,ndata)
  qdferr(2,i,*) =  reform(cos(laterr(i,*))*sin(azierr(i,*)),1,1,ndata)
  qdferr(1,i,*) =  reform(cos(laterr(i,*))*cos(azierr(i,*)),1,1,ndata)
  
  qdferr(*,i,*) = reform(q_rot(q_conj(qrot),reform(qdferr(*,i,*),4,ndata)),4,1,ndata)
endfor

;thn = acos(qdfn(3,*))
;phn = atan(qdfn(2,*),qdfn(1,*))
;ths = acos(qdfs(3,*))
;phs = atan(qdfs(2,*),qdfs(1,*))
;the = acos(qdfe(3,*))
;phe = atan(qdfe(2,*),qdfe(1,*))
;thw = acos(qdfw(3,*))
;phw = atan(qdfw(2,*),qdfw(1,*))

end