Viewing contents of file '../idllib/iuedac/iuelib/pro/barker.pro'
;*************************************************************************
;+
;*NAME:
;
; BARKER
;
;*CLASS:
;
; Spectral Data Reduction
;
;*CATEGORY:
;
;*PURPOSE:
;
; To derive the ripple parameters K and alpha for orders
; m1 to m2 for IUE high dispersion data. The algorithm
; used is based on the Barker ripple correction algorithm
; A.J. 89, 899 (1984), which relies upon minimizing the
; flux differences in order overlap region.
;
;*CALLING SEQUENCE:
;
; BARKER,IMAGET,AIR,VRAD,M1,M2,A,K1,K,DK,ER,NCAM,ISN
;
;*PARAMETERS:
;
; IMAGET (REQ) (I)
; IUE G.O. file name
;
; AIR (REQ) (I)
; flag for supressing conversion to vacuum wavelengths
; AIR=1 will keep air wavelengths
; AIR=0 will go to vacuum wavelengths (w>2000. only)
;
; VRAD (REQ) (I)
; Radial velocity specified to correct SIPS wavelength
; scale. This parameter should be non-zero only for
; IUE G.O. files processed under the "old" software,
; for which no heliocentric velocity correction was
; applied.
;
; M1 (REQ) (I)
; starting order number
;
; M2 (REQ) (I)
; final order number N.B. code assumes M2>M1
;
; A (REQ) (O)
; Ahmad ripple fudge factor for IUE
;
; K (REQ) (I)
; ripple parameters for orders M1 to M2
;
; DK (REQ) (O)
; increment in K (in units of K)
;
; ER (REQ) (O)
;
; NCAM (REQ) (I)
; camera number
;
; ISN (REQ) (I)
; image sequence number
;
;*EXAMPLES:
;
;*SYSTEM VARIABLES USED:
;
; !PI
;
;*INTERACTIVE INPUT:
;
;*SUBROUTINES CALLED:
;
; PARCHECK
; FILETYPE
; RIPGET
;
;*FILES USED:
;
;*SIDE EFFECTS:
;
;*RESTRICTIONS:
;
; The SIPS implementation of the ripple correction assumes
; (see SIPS V. 2. 6.4.2) that the heliocentric velocity correction
; has been done before the ripple correction. At this time, the
; wavelength scale is still in vacuum wavelengths. RDAF applications
; of the ripple correction work from the G.O. files which have
; all wavelengths longward of 2000 A converted to AIR wavelengths.
;
;*NOTES:
;
; tested with IDL Version 2.1.0 (sunos sparc) 23 Jul 91
; tested with IDL Version 2.1.0 (ultrix mispel) N/A
; tested with IDL Version 2.1.0 (vms vax) 23 Jul 91
;
;*PROCEDURE:
;
; The Ahmad ripple parameter A is held constant for each
; camera, and the K parameter is allowed to vary.
;
;*I_HELP nn:
;
;*MODIFICATION HISTORY:
;
; Programmer: T.B. Ake (1983)
; Version 0 (IUE Calibration) T.B. Ake 1983
; Version 1 CAG transfer to VAX, add documentation
; add error checking.
; Version 2
; 1-09-91 PJL modified to sun/unix; included removing !ERR
; 7-23-91 PJL cleaned up; added PARCHECK; tested on SUN and VAX;
; updated prolog
; 9-30-91 PJL corrcted typo effecting lwp's a
;
;-
;************************************************************************
pro barker,imaget,air,vrad,m1,m2,a,k1,k,dk,er,ncam,isn
;
; parse file name, verify that data file is on disk, and get camera number
;
npar=n_params(0)
if npar eq 0 then begin
print,'BARKER,IMAGET,AIR,VRAD,M0,M1,A,K,DK,ER,NCAM,ISN'
return
endif ; npar
parcheck,npar,12,'BARKER'
filetype,imaget,disk,path,name,ext,vers,ncam,type,isn
;
; only allow barker ripple correction on IUE echelle data
;
if (strlowcase(type) ne 'h') then begin
print,'File ',imaget, ' is not a high dispersion extracted spectral file.'
print,'ACTION: Returning to calling procedure.'
return
endif ; type
;
; define the ripple parameters k for each of the orders
; for which ripple constants are desired
;
k=fltarr(m2+1-m1)
dk=k
er=k
;
; get default values for k0 and a based on the camera number
;
case 1 of
(ncam eq 1) : begin ; lwp
k0=231150. ;start with same set of constants
;input parameter
a=0.896 ;for both the lwp and the lwr (guess)
minm=75
maxm=124
end ; ncam eq 1
(ncam eq 2) : begin ; lwr
k0=231150.
a=0.896
minm=77
maxm=124
end ; ncam eq 2
(ncam eq 3) : begin ; swp
k0=137750.
a=0.856
minm=62
maxm=124
end ; ncam eq 3
else: begin ; not a useful camera
print,' Invalid camera number requested'
print,' ACTION: Returning to calling procedure'
return
end ; else
endcase ; 1
;
; check user-specified order numbers for valid range for each camera
;
if (m2 gt maxm) or (m1 lt minm) then begin
print, 'Highest or lowest order number is out of range for camera ',ncam
print, 'ACTION: Returning to calling procedure'
return
endif ; m2
;
; now that the error checking is out of the way, open the data file for reading
;
openr,unit,/get_lun,disk + path + name + '.dat' + vers, error=err
;
; read in the initial order
;
ripget,unit,m1-1,air,vrad,k0,a,w,f,x
;
; identify wavelength range to be considered
;
minw=where(((x ge 0.20) and (w lt k0/(m1-1))),count)
if count lt 1 then minw=w(0) else minw=w(minw(0))
second=0 & avg=fltarr(4) & tmp=avg
;
; go from order m1 to order m2 (e.g. begin at longer wavelengths)
;
for i=m1,m2 do begin
again:
m=i+second
ripget,unit,m,air,vrad,k0,a,w1,f1,x
;
; identify the region of wavelength overlap
;
maxw=where(((x lt 0.20) and (w1 gt k0/m)),count)
if count lt 1 then maxw=max(w1) else maxw=w1(maxw(0))
if minw gt maxw then begin
print,' order',m-1,',',m,' no order overlap'
goto,next
endif ; minw
;
; get the indices for the wavelengths in the order overlap region
;
inx=where(((w ge minw) and (w le maxw)),count)
;
; use system variable to get number of points in overlap wavelength
; range for the first order
;
np = count
inx1=where(((w1 ge minw) and (w1 le maxw)),count)
;
; find number of points in the overlap wavelength range for the second
; order using system variable.
;
nm = count
k1=k0 & dk1=1. & iter=0
pia=!pi*a ; (treated as a constant for this program)
;
; restrict number of iterations to 50 (maximum)
;
while ((abs(dk1) gt .5) and (iter lt 50)) do begin
; print,'iteration number ',iter
;
; estimate the ripple corrected flux for the first order
;
xm=pia*(m-1)*(1-k1/(m-1)/w(inx))
r=sin(xm)/xm
fp=f(inx)/r/r
;
; estimate the ripple corrected flux for the second order
;
xm1=pia*m*(1-k1/m/w1(inx1))
r1=sin(xm1)/xm1
fm=f1(inx1)/r1/r1
;
; erf is effectively the ratio of the ripple corrected fluxes,
; normalized by the number of points for each order in the
; wavelength overlap region.
;
erf=total(fp)/np/total(fm)*nm
;; erf=total(2*abs(fp-fm))/(abs(fp)+abs(fm))
drm=-2/pia/((m-1)*w(inx)-k1)*(cos(xm)-pia*r)
drm1=-2/pia/(m*w1(inx1)-k1)*(cos(xm1)-pia*r1)
;
; compute the increment in k1 based on barker criterion
;
dk1=(erf-1)/(total(drm)/np-erf*total(drm1)/nm)
k1=k1+dk1 & iter=iter+1
endwhile ; dk1
print,'order =', m,' iteration number ',iter
tmp(0)=1 & tmp(1)=k1 & tmp(2)=(k1-k0)*(k1-k0) & tmp(3)=erf
avg=avg+tmp
next:
w=w1 & f=f1
minw=where(((x ge 0.20) and (w lt k0/m)),count)
if count lt 1 then minw=w(0) else minw=w(minw(0))
if second eq 0 then begin
second=1
goto,again
endif ; second
j=i-m1
;
; load the ripple data into output variables
;
nord=avg(0)>1
k(j)=avg(1)/nord
if nord gt 1 then dk(j)=sqrt(avg(2)-(avg(1)-nord*k0)*(avg(1)-nord*k0)/nord)
er(j)=avg(3)/nord
avg=tmp
end ; i
free_lun,unit
return
end ; barker