Viewing contents of file '../idllib/contrib/icur/addheliocor.pro'
addheliocor.pro;2 000640 000451 000017 00000000000 06615503463 015063 0 ustar 00fwalter users 000000 000000 addline.pro 000640 000451 000017 00000003257 06615503076 014104 0 ustar 00fwalter users 000002 247506 ;*************************************************************
PRO ADDLINE,WAVE,FLUX,lam,sig,tf,helpme=helpme ; ADD LINE TO FLUX VECTOR
COMMON COM1,H
common comxy,xcur,ycur,zerr,resetscale,lu3
if keyword_set(helpme) then begin
print,' '
print,'* ADDLINE,WAVE,FLUX,lam,sig,tf '
print,'* Add a single line to the flux vector'
print,'* WAVE, FLUX: wavelength and flux vectors. FLUX is modified'
print,'* LAM,TF: wavelength and integrated flux of line'
print,'* SIG: sigma in bins or -Angstroms (def=2)
print,'* '
print,'* The LAM, SIG, and TF values are requested if not passed.'
print,' '
return
endif
S=N_ELEMENTS(WAVE)
DISP=(WAVE(50)-WAVE(0))/50.
if n_params() ge 5 then pass=1 else pass=0
if pass then begin
if (lam le 0) or (tf le 0.) then pass=0
endif
if not pass then begin
LAM=0. & SIG=1. & TF=1.E-15
z=' ENTER WAVELENGTH, SIGMA (BINS), AND INTEGRATED FLUX: '
;print,z
READ,LAM,SIG,TF,prompt=z
endif
IF TF EQ 0. THEN RETURN
if lam le 0. then return
if sig eq 0. then sig=2. ;def=2 bins
if sig lt 0. then sig=abs(sig/disp) ;passed in Angstroms
AMP=TF/(SIG*DISP*2.5066)
G=FLTARR(S)
i0=xindex(wave,lam) ;TABINV,WAVE,LAM,I0
WID=FIX(SIG)*4
I1=FIX(I0)-WID
IF I1 LT 0 THEN I1=0
I2=FIX(I0)+WID
IF I2 GE S THEN I2=S-1
FOR I=I1,I2 DO BEGIN
Z=(FLOAT(I)-I0)/SIG
Z=Z*Z
G(I)=G(I)+AMP*EXP(-Z/2.)
ENDFOR
FLUX=FLUX+G
if n_elements(h) lt 65 then return
H(60)=999
H(61)=H(61)+1
J=62+(H(61)-1)*3
H(J)=FIX(LAM)
H(J+1)=FIX(1000.*(LAM-FIX(LAM)))
H(J+2)=FIX(ALOG10(ABS(TF))*100.)
if n_elements(lu3) le 0 then begin
printf,lu3,' 6'
printf,lu3,lam,sig,TF,' ADDLIN: Lambda,Sigma,Flux'
endif
RETURN
END
addline.pro;2 000640 000451 000017 00000000000 06615503076 016361 1addline.pro ustar 00fwalter users 000002 247506 addred.pro 000640 000451 000017 00000011255 06615503077 013713 0 ustar 00fwalter users 000002 253310 ;*****************************************************************************
PRO ADDRED,IM,WAVE,FLUX,EBMV,trans ; CORRECT FOR INTERSTELLAR EXTINCTION
; TAKEN FROM [210,021]NUNRED.PRO 4/5/83
; IM=0,-1 FOR UNRED
COMMON COM1,H
COMMON VARS,VAR1,VAR2,VAR3,VAR4,VAR5
COMMON ICDISK,ICURDISK,ICURDATA,ISMDAT
common comxy,xcur,ycur,zerr,rsc,lu3
if n_elements(var3) eq 0 then var3=(2*256) ;REDDENING TABLE
nhd=n_elements(h)
TABLE=FIX(VAR3 AND 7*256)/256
TABLE=FIX(TABLE)
IF (TABLE LT 0) OR (TABLE GT 5) THEN TABLE=2
IF N_PARAMS(0) LE 3 THEN begin
IF N_ELEMENTS(H) GT 92 THEN BEGIN
if H(92) ne 0 then cext=string(float(H(92))/1000.,'(F6.2)') else cext=''
z=' Current E(B-V)='+cext
if strlen(cext) gt 0 then print,z
endif
READ,' Enter E(B-V); - to redden: ',EBMV
endif
EBMV=-FLOAT(EBMV)
trans=wave*0.+1.
IF EBMV EQ 0. THEN begin
if nhd gt 92 then begin
h(91)=0
h(92)=0
endif
RETURN
endif
IF ABS(EBMV LT 100.) THEN NH=EBMV*5.9E21 ELSE NH=EBMV
Q=912.
SI=N_ELEMENTS(WAVE)-1
SIG=WAVE*0.
LMAX=WAVE(SI)
APB=(LMAX-WAVE(0))/FLOAT(SI)
IMIN=(((Q-WAVE(0))/APB > 0) < SI)
IF LMAX GT Q THEN BEGIN ; USE A TABULATED REDDENING LAW
;
; if user is to enter the name of the table file
GET_LUN,LUN
TABFILE=''
IF TABLE EQ 0 THEN REPEAT BEGIN
PRINT,'Please enter the name of the file containing the extinction curve'
PRINT,'you want to use. The wavelength information must be in record one'
PRINT,'and the flux in record 2.'
READ,'What is the name of the file ?',TABFILE
OPENR,LUN,TABFILE
END UNTIL !ERR GT 0
;
IF TABLE EQ 1 THEN TABFILE=ismdat+'SAVMAT.DAT' ; Savage and Mathis
IF TABLE EQ 2 THEN TABFILE=ismdat+'SEATON.DAT' ; Seaton
IF TABLE EQ 3 THEN TABFILE=ismdat+'NANDY.DAT' ; Nandy
IF TABLE EQ 4 THEN TABFILE=ismdat+'ORI.DAT' ; Bohlin & Savage
IF TABLE EQ 5 THEN TABFILE=ismdat+'SMC.DAT' ; Sm Magellenic Clouds
OPENR,LUN,TABFILE
tlun=fstat(lun) & zlen=tlun.rec_len/4
Z=ASSOC(LUN,FLTARR(zlen)) ; Z is always associated variable
XTAB=Z(0) ; S & M wavelength in record 0
YTAB=Z(1) ; S & M flux in record 1
CLOSE,LUN & FREE_LUN,LUN
I=WHERE(XTAB GT 0)
XTAB=FLOAT(XTAB(I)) ; extract valid points and make sure they
YTAB=FLOAT(YTAB(I)) ; are floating point format
F=FLOAT(FLUX) ; make sure input wave and flux are
W=FLOAT(WAVE) ; floating point format
TAB='??!! SEE DATA AID !!??'
IF TABLE EQ 1 THEN TAB=' (Savage and Mathis) '
IF TABLE EQ 2 THEN TAB=' (Seaton) '
IF TABLE EQ 3 THEN TAB=' (Nandy) '
IF TABLE EQ 4 THEN TAB=' (Orion) '
IF TABLE EQ 5 THEN TAB=' (SMC) '
IF TABLE EQ 0 THEN TAB=TABFILE
;PRINT,'The flux is being unreddened.',TAB,'E(B-V) = ',EBMV
X=INTERPOL(YTAB,XTAB,W) ; interp table to current wavelngth scale
SIG(IMIN)=1.65E-22*(X(IMIN:*)+3.1)
LMAX=Q
ENDIF
IF WAVE(0) LT Q THEN BEGIN
FOR I=0,IMIN DO BEGIN
LAM=WAVE(I)
IF (LAM LE 912.) AND (LAM GT 504.) THEN BEGIN ; H-He
SLO=2.58 & INT=-24.86 & GOTO,SIGM
ENDIF
IF (LAM LE 504.) AND (LAM GT 43.648) THEN BEGIN ; He-C
SLO=2.72 & INT=-25.05 & GOTO,SIGM
ENDIF
IF (LAM LE 43.648) AND (LAM GT 30.99) THEN BEGIN ; C-N
SLO=2.87 & INT=-25.22 & GOTO,SIGM
ENDIF
IF (LAM LE 30.99) AND (LAM GT 23.301) THEN BEGIN ; N-O
SLO=3.06 & INT=-25.48 & GOTO,SIGM
ENDIF
IF (LAM LE 23.301) AND (LAM GT 14.3018) THEN BEGIN ; N-Ne
SLO=2.54 & INT=-24.52 & GOTO,SIGM
ENDIF
IF (LAM LE 14.3018) AND (LAM GT 9.5122) THEN BEGIN ; Ne=Mg
SLO=2.52 & INT=-24.46 & GOTO,SIGM
ENDIF
IF (LAM LE 9.5122) AND (LAM GT 6.738) THEN BEGIN ; Mg-Si
SLO=1.91 & INT=-23.80 & GOTO,SIGM
ENDIF
IF (LAM LE 6.738) AND (LAM GT 5.019) THEN BEGIN ; Si-S
SLO=3.73 & INT=-25.27 & GOTO,SIGM
ENDIF
IF (LAM LE 5.019) AND (LAM GT 3.871) THEN BEGIN ; S-A
SLO=2.96 & INT=-24.68 & GOTO,SIGM
ENDIF
IF LAM LE 3.871 THEN BEGIN ;FINALLY SHORTWARD OF ARGON K EDGE
SLO=2.99 & INT=-24.68 & GOTO,SIGM
ENDIF
SIGM:SIG(I)=10.^(((SLO)*ALOG10(LAM))+INT)
ENDFOR
ENDIF
T=(SIG*NH<80.)
TRANS=EXP(-T)
FLUX=Flux*TRANS ; F*10.^(.4*EBMV*(X+3.1)) ; derive unreddened flux
if (im ne -1) and (nhd gt 0) and (n_elements(lu3) gt 0) then begin
H(91)=TABLE ; unreddening flag
H(92)=H(92)+FIX(1000*EBMV) ; total E(b-v)*1000
IF H(92) EQ 0 THEN H(91)=0 ;NO REDDENING
printf,lu3,'-5'
printf,lu3,'Data Unreddened: Table=',TABLE,' ',TAB,'; E(B-V)=',EBMV
endif
RETURN
END
addred.pro;2 000640 000451 000017 00000000000 06615503077 016016 1addred.pro ustar 00fwalter users 000002 253310 addsp.pro 000640 000451 000017 00000017611 06615503077 013565 0 ustar 00fwalter users 000002 260510 ;**************************************************************************
pro addsp,file,h0,w,f0,e0,prlog=prlog,debug=debug
common trim,ntr,ntr1,range,kr,irg,r2
;
;
swlw=1255. & swuw=1265. ;Si
lwlw=2794. & lwuw=2797. ;MgII k
print,' ADDSP : procedure to coadd spectra'
print,' set debug=-1 to turn off cross correlation and add with no bin shift'
print,' set debug=2 to run fully interactively and query about all shifts'
print,' '
shftmode=1
if n_params(0) eq 0 then file=''
if file eq '' then read,' enter name of data file: ',file
!quiet=1
;
if not ffile(file) then begin
print,' data file not found: returning
return
endif
;
get_lun,lu
print,' these are the current data file contents:'
ldat,file
read,' enter record number for storage, -1 to append, -9 to skip: ',savrec
if (savrec lt -1) and (n_params(0) lt 4) then begin
print,' WARNING: no data being saved or passed to MAIN'
print," To pass data back, call ADDSP,'',H,W,F,E"
return
endif
savrec=fix(savrec)
title=''
read,' enter title for header: ',title
;
free_lun,lu
print,' what records do you wish? (or first, -#) -1 to end'
print,' Note that the first record is the template, and should be properly exposed'
read,recs
if recs eq -1 then return
a=0
while a ge 0 do begin
read,a
recs=[recs,a]
endwhile
nr=n_elements(recs)
if (nr eq 2) and (recs(1) lt 0) then begin
nr=abs(recs(1))+1
recs=recs(0)+indgen(nr)
endif
nr=fix(nr)-1
recs=fix(recs(0:nr-1)) ;cut last value
print,nr,' Records chosen:',recs
if nr le 1 then stop
if not keyword_set(debug) then debug=0
if debug eq -1 then iccr=0 else iccr=1
if debug ne -1 then begin
read,' if you wish to cross correlate the spectra, enter 1',iccr
print,' You have 2 choices: automatic cross correlations (0) or'
print,' manual shifting (default) (1)'
read,' Enter your choice: ',iccr2
if iccr2 ne 0 then iccr2=1
print,' '
eta=(nr-1)*7+3
; if iccr eq 1 then print,'*** Estimated time to completion is',eta,' minutes ***'
endif
;
get_lun,lu
openw,lu,'addsp.log'
print,' '
printf,lu,' ADDSP processing log'
print,' processing starting at ',systime(0)
printf,lu,' ADDSP processing starting at ',systime(0)
printf,lu,nr,' Records chosen:',recs
printf,lu,' Title:',title
if iccr eq 1 then z=' ' else z=' not '
printf,lu,' Spectra will',z,'be cross correlated'
print,' '
gdat,file,h0,w0,f0,e0,recs(0)
ncam=h0(3)
image=h0(4)
if (ncam le 4) and (image lt 0) then image=image+65636L
if ncam le 4 then printf,lu,' this camera is number ',ncam,' image',image $
else printf,' Data is ',strtrim(byte(h(100:159)>32),2)
addspwl,h0,w0,lam1,lam2,dl,krange
printf,lu,' Wavelengths limits are ',lam1,' to',lam2,' ; increment=',dl
printf,lu,' '
w=lam1+dl*findgen((lam2-lam1)/dl)
e0=fix(e0)
finter,w,w0,f0,e0
if (ncam le 4) and (n_elements(w) gt 4000) then kgap,h0,w,wgap,ngap,e0
fref=f0
eref=e0 ;
if (ncam le 4) and (n_elements(w) gt 4000) then begin
ksplice,h0,w,eref
; fftsm,fref,1,0.6
endif
f0=f0*(h0(5)>1) ;multiply by observing time
time=long(h0(5))>1L
h0(180)=h0(4)
k=wherebad(e0,1) ;bad data
f0(k)=0.
e0=long(e0)*0+(long(h0(5))>1L)
e0(k)=0
;
nrecs=nr-1
if nr eq 2 then nit=-1 else nit=0
for irecs=1,nrecs do begin
gdat,file,h1,w1,f1,e1,recs(irecs)
im2=h1(4)
if (h1(3) le 4) and (im2 lt 0) then im2=im2+65636L
if h1(3) ne ncam then begin
print,' Warning: camera =',h1(3),' not ',ncam
printf,lu,' Warning: camera =',h1(3),' not ',ncam
if ((ncam eq 1) and (h1(3) eq 2)) or ((ncam eq 2) and (h1(3) eq 1)) then $
goto, camok
stop
irecs=irecs+1
nr=nr-1
goto,done
endif
camok:
if ncam le 4 then printf,lu,' Processing camera number',h1(3),', image',im2 $
else printf,lu,' Processing record ',strtrim(byte(h(100:159)>32b),2)
e1=fix(e1)
finter,w,w1,f1,e1
if (ncam le 4) and (n_elements(w) gt 4000) then begin
kgap,h1,w,wg1,ng,e1
for i=0,ng-1 do begin
if i lt ngap then if abs(wg1(i,0)-wgap(i,0)) lt 3. then begin
wgap(i,0)=(wgap(i,0)>wg1(i,0))
wgap(i,1)=(wgap(i,1)<wg1(i,1))
endif ;reset gaps
endfor ;gap loop
endif ;idat = 7 gap loop
e2=e1
if (ncam le 4) and (n_elements(w) gt 4000) then begin
ksplice,h1,w,e2
; fftsm,f1,1,0.6
endif
;
if iccr eq 1 then begin
if debug ne 0 then print,' beginning SPSHFT'
nit=nit+1
s=spshft(fref,f1,eref,e2,nit)
if iccr2 eq 0 then begin ; autoshift
wk=where(w gt lw) & kw1=(wk(0)-150)>0
wk=where(w gt uw) & kw2=(wk(0)+150)<(n_elements(w)-1)
print,' bin limits=',kw1,kw2
print,systime(0)
s=spshft(fref(kw1:kw2),f1(kw1:kw2),eref(kw1:kw2),e2(kw1:kw2),nit)
print,' SPSHFT done, s=',s,systime(0)
endif else begin ;manual shift
s=shfth2(w,fref,f1,lw,uw)
print,' SHFTH2 done, s=',s
;s=0
endelse ;iccr2
endif else s=0. ;iccr
;
igo=0
case 1 of
s eq -1000: begin
printf,lu,' WARNING: file number',irecs+1,': peak too close to edge'
end
s lt -8000: begin
printf,lu,' file number',irecs+1,' WARNING: NO SIGNIFICANT PEAK FOUND'
s=s+9000
igo=1
end
abs(s) gt krange: begin
z=string(s,'(I2)')
printf,lu,' file number',irecs+1,' WARNING: PEAK @ ',z,' OUTSIDE BOUNDS - RESETTING'
igo=1
end
else:
endcase
if (igo eq 1) or (debug eq 2) then begin
if debug gt 0 then begin ;check manually
print,string(7b),' shift=',s
igo=0.0
read,' is this shift OK? enter 99 to accept,-99 to omit, other to set: ',igo
case 1 of
igo ge 98.9:
igo le -98.9: s=-1000
else: s=igo
endcase
endif else s=-1000 ;ignore data if running on auto
endif ;igo
if s eq -1000 then begin ;skip this record
printf,lu,' file number',irecs+1,' not added'
nr=nr-1 & goto,done
endif
f1=f1*(h1(5)>1)
time=time+(long(h1(5))>1L)
h0(180+irecs)=h1(4)
k=wherebad(e1,1) ;bad data
f1(k)=0.
e1=long(e1)*0+(long(h1(5))>1L)
e1(k)=0
if s ne 0 then begin ;apply shift
if (shftmode eq 0) or (abs(s-fix(s)) lt 0.01) then begin
s=fix(s)
f1=shift(f1,s)
e1=shift(e1,s)
case 1 of
s gt 0: begin
f1(n_elements(f1)-1-s:*)=0.
e1(n_elements(f1)-1-s:*)=0
end
s lt 0: begin
f1(0:abs(s)-1)=0.
e1(0:abs(s)-1)=0
end
else:
endcase
endif else begin
dl=w(1)-w(0)
shft=s*dl
f1=interpol(f1,w,w+shft)
e1=interpol(e1,w,w+shft)
endelse
endif
f0=f0+f1
e0=e0+e1
print,' file number',irecs+1,' added, shift=',s,' accumulated time=',time
printf,lu,' file number',irecs+1,' added, shift=',s,' accumulated time=',time
done:
endfor ;irecs
f0=f0/(e0>1) ;flux/sec
;
if (ncam le 4) and (n_elements(w) gt 4000) then cutgaps,h0,w,f0,e0,ngap,wgap(*,0),wgap(*,1)
;
h0(4)=0 ;dummy record number
if time gt 32767L then h0(5)=-fix(time/60) else h0(5)=time
h0(9)=nr ;records added
title=title+' sum of '+strtrim(nr,2)+' spectra'
title=byte(title)
h0(100)=title
;
print,h0(20:23)
e0=byte(127.*float(e0)/max(e0))
if savrec ne -1 then begin
kdat,file,h0,w,f0,e0,savrec
printf,lu,' data saved to ',file,' record #',savrec
endif
!x.title='!6 Angstroms'
!y.title=ytit(0)
!p.title=strtrim(byte(title),2)
!p.font=-1
!p.charsize=7./5. ;!fancy=3
print,' '
printf,lu,' '
print,' Done at ',systime(0)
printf,lu,' ADDSP done at ',systime(0)
close,lu
free_lun,lu
if !version.os eq 'vms' then z='delete/noconfirm spshft.tmp;*' $
else z='rm spshft.tmp'
spawn,z
if keyword_set(prlog) then spawn_print,'addsp.log' else print,' Log in addsp.log'
print,' '
return
end
addsp.pro;2 000640 000451 000017 00000000000 06615503077 015536 1addsp.pro ustar 00fwalter users 000002 260510 addspec.pro 000640 000451 000017 00000002374 06615503077 014116 0 ustar 00fwalter users 000002 247576 ;**************************************************************************
pro addspec,w,f,w1,f1,e0,w0,f0
w0=w
f0=f
nw0=n_elements(w0)
nw1=n_elements(w1)
dw=(w0(nw0-1)-w0(0))/(nw0-1)
ov1=where((w1 gt w0(0)) and (w1 lt w0(nw0-1))) ;overlap region
ov0=where((w0 gt w1(0)) and (w0 lt w1(nw1-1))) ;overlap region
f1=f1*mean(f(ov0))/mean(f1(ov1)) ;set scaling
if n_elements(e0) eq nw0 then f0=f*e0 else e0=intarr(nw0)+1
;
if w1(0)+dw lt w0(0) then begin ;second vector starts first
np1=(w0(0)-w1(0))/dw
w0=[w0(0)-np1*dw+dw*findgen(np1),w0]
endif else np1=0
nw0=n_elements(w0)
;
if w1(nw1-1)-dw gt w0(nw0-1) then begin ;second vector ends last
np2=(w1(nw1-1)-w0(nw0-1))/dw
w0=[w0,w0(nw0-1)+dw+dw*findgen(np2)]
endif else np2=0
;
np1=fix(np1+0.5)
np2=fix(np2+0.5)
np=n_elements(w0)
e=intarr(np)
f2=fltarr(np)
f2(np1)=f0
e(np1)=e0
;
f3=interpol(f1,w1,w0)
k=where(w0 le w1(0)+dw)
if k(0) ge 0 then begin
f3(k)=0.
km=max(k)+1
endif else km=0
k=where(w0 ge w1(nw1-1)-dw)
kgood=k(0)
if kgood gt -1 then f3(k)=0.
if kgood eq -1 then kgood=n_elements(f3)-1
np=kgood-km ;+1
k=km+indgen(np)
e(k)=e(k)+1
f0=f2+f3
k=where(e gt 0)
e=e(k)
f0=f0(k)
w0=w0(k)
f0=f0/e
e0=e
if n_params(0) lt 7 then begin
w=w0
f=f0
endif
return
end
addspec.pro;2 000640 000451 000017 00000000000 06615503077 016377 1addspec.pro ustar 00fwalter users 000002 247576 addspwl.pro 000640 000451 000017 00000002544 06615503077 014130 0 ustar 00fwalter users 000002 260322 ;****************************************************************
pro addspwl,h0,w0,lam1,lam2,dl,krange ;get wavelength limits
idat=h0(0)
ncam=h0(3)
case 1 of
idat eq 0: begin ;IUE low dispersion
krange=7
case 1 of
ncam eq 3: begin ;SWP
lam1=1100.
lam2=1980.
dl=1.0
end
(ncam eq 1) or (ncam eq 2): begin ;LWP/R
lam1=1900.
lam2=3300.
dl=1.5
end
else: begin
print,' invalid camera number =',ncam
stop
end
endcase
end
idat eq 1: begin
krange=11
lam1=min(w0)
lam2=max(w0)
dl=float(h(22))+float(h(23))/1000.
end
idat eq 7: begin ;IUE high dispersion, long format
krange=11
case 1 of
ncam eq 3: begin ;SWP
lam1=1150.
lam2=2068.5
end
(ncam eq 1) or (ncam eq 2): begin ;LWP/R
lam1=2000.
lam2=3200.
end
else: begin
print,' invalid camera number =',ncam
stop
end
endcase
dl=(w0(5000)-w0(0))/5000.
end
else: begin
print,' ADDSPWL: parameters for IDAT=',idat,' undefined'
krange=11
lam1=min(w0)
lam2=max(w0)
dl=w0(1)-w0(0)
end
endcase
return
end
addspwl.pro;2 000640 000451 000017 00000000000 06615503077 016445 1addspwl.pro ustar 00fwalter users 000002 260322 amod.pro 000640 000451 000017 00000003067 06615503100 013405 0 ustar 00fwalter users 000002 252265 ;*********************************************************
PRO AMOD,A,ISHAPE
;*********************************************************
;
PARAM=-1 & FACTOR=1. & OPERATION='*'
REPEAT BEGIN
IF (PARAM EQ -1) THEN BEGIN
PRINT,'Here is the A vector, as it now stands:'
PRINT,' '
PRINT,' INTENSITY POSITION WIDTH'
FORM='("LINE #",I2,11X,E11.4,F10.2,F9.2)'
FOR I=0,ISHAPE-1 DO PRINT,format=form,I,A(3*I),A(3*I+1),A(3*I+2)
PRINT,' '
PRINT,'YOU will be prompted for the PARAMeter to fudge and'
PRINT,'then the fudge FACTOR and OPERATION (+,-,*,/,=)'
PRINT,' '
PRINT,'Enter -9 to RETURN; ENTER -1 to display current A vector'
PRINT,' '
PRINT,' '
ENDIF ELSE BEGIN
IF (PARAM GE 3*ISHAPE) OR (PARAM LT -1) THEN BEGIN
PRINT,'NO SUCH PARAMETER, TRY AGAIN'
PARAM=-1
ENDIF ELSE BEGIN
PRINT,' Original PARAMeter (',PARAM,') = ',A(PARAM)
READ,'fudge FACTOR? ',FACTOR
READ,' OPERATION? ',OPERATION
OPERATION=STRTRIM(OPERATION)
IF OPERATION EQ '*' THEN A(PARAM)=A(PARAM)*FACTOR
IF OPERATION EQ '/' THEN A(PARAM)=A(PARAM)/FACTOR
IF OPERATION EQ '+' THEN A(PARAM)=A(PARAM)+FACTOR
IF OPERATION EQ '-' THEN A(PARAM)=A(PARAM)-FACTOR
IF OPERATION EQ '=' THEN A(PARAM)=FACTOR
PRINT,' FUDGED PARAMeter (',PARAM,') = ',A(PARAM)
ENDELSE
ENDELSE
READ,'PARAMETER? ',PARAM
ENDREP UNTIL (PARAM EQ -9)
;
PRINT,'continue entering parameters to fix... (-9 when finished; -10 to stop)'
RETURN
END
amod.pro;2 000640 000451 000017 00000000000 06615503100 015203 1amod.pro ustar 00fwalter users 000002 252265 angstrom.pro 000640 000451 000017 00000000065 06615503100 014314 0 ustar 00fwalter users 000002 252267 function angstrom,dum
return,'!6!sA!r!u!9 %!6!n'
end
angstrom.pro;2 000640 000451 000017 00000000000 06615503100 017031 1angstrom.pro ustar 00fwalter users 000002 252267 appendicd.pro 000640 000451 000017 00000000471 06615503100 014405 0 ustar 00fwalter users 000002 253315 ;****************************************************************************
pro appendicd,file1,file0
icdfile,file1,nr
for i=0,nr do begin
gdat,file1,h,w,f,e,i
if n_elements(h) eq 1 then goto,endfile
kdat,file0,h,w,f,e,-1
endfor
endfile:
print,i,' files copied from ',file1,' to ', file0
return
end
appendicd.pro;2 000640 000451 000017 00000000000 06615503100 017216 1appendicd.pro ustar 00fwalter users 000002 253315 avbad.pro 000640 000451 000017 00000002763 06615503100 013536 0 ustar 00fwalter users 000002 252700 ;***********************************************************************
function avbad,f1,e,kk,bad=bad ;expunge bad data points
;
f=f1
; first, get rid of single points
f=avspt(f,e,kz,bad=bad)
if n_elements(kz) eq 0 then kk=-1 else kk=kz ;indices of single pts
;
np=n_elements(f)-1
if n_elements(bad) gt 0 then k=bad else case 1 of
(min(e) eq 0) and (max(e) lt 10): k=where(e gt 1)
(min(e) lt 0) : k=where(e lt -200)
else : k=where(e eq 0)
endcase
nk=n_elements(k) ;number of bad points
nkk=n_elements(kz) ;number of bad single points
nsp=nk-nkk ;number of other bad points
if nk le 0 then return,f
if nsp le 0 then return,f
;
if nkk gt 0 then begin ;there are single bad points
for i=0,nkk-1 do k(where(k eq kk(i)))=-1
k=k(where(k gt -1))
nk=n_elements(k)
if nk ne nsp then begin
print,' AVBAD: Warning: nk=',nk,' nk-nkk=',nsp
endif
endif
;
k1=k(1:*)-k ;difference
k2=where(k1 gt 1,nk2)
if nk2 eq 0 then k2=[0,n_elements(k)] else k2=[0,k2+1,n_elements(k)]
nk2=nk2+1 ;number of bad intervals
print,nk2,' bad intervals being interpolated'
for i=0,nk2-1 do begin
i1=k(k2(i))-1 & i2=k(k2(i+1)-1)+1
case 1 of
i1 eq -1: f(0:i2-1)=f(i2)
i2 ge np: f(i1+1:*)=f(i1)
else: begin
np1=i2-i1 ;number of points
sl=(f(i2)-f(i1))/np1
x=f(i1)+findgen(np1)*sl
f(i1)=x
end
endcase
endfor
;
return,f
end
avbad.pro;2 000640 000451 000017 00000000000 06615503100 015467 1avbad.pro ustar 00fwalter users 000002 252700 averagespec.pro 000640 000451 000017 00000004062 06615503101 014745 0 ustar 00fwalter users 000002 251444 ;*******************************************
function AVerageSPEC,FILE,i1,i2,w,save=save,helpme=helpme,noweight=noweight
if n_params(0) lt 2 then helpme=1
if keyword_set(helpme) then begin
print,' '
print,'* AVERAGESPEC - average spectra in .ICD file '
print,'* calling sequence: f=averagespec(file,i1,i2,w) or f=averagespec(file,rec,w)'
print,'* default uses time weighting'
print,'* FILE: name of .ICD file'
print,'* I1,I2: first and last records to be averaged'
print,'* REC: list of records to be averaged'
print,'* W: output wavelength vector (optional)'
print,'* '
print,'* KEYWORDS'
print,'* NOWEIGHT: if set, do straight average, default weights by times'
print,'* SAVE: if set, append average spectrum to end of input .ICD file'
print,' '
return,0
end
;
ext=get_ext(file)
if strlen(ext) eq 0 then ext='.icd'
if not ffile(file+ext) then begin
ic=getenv('icurdata')
if ffile(ic+file+ext) then file=ic+file else begin
print,' file ',file,' not found - returning'
return,0
endelse
endif
if n_elements(i1) gt 1 then irec=1 else irec=0
if irec then k=i1 else k=i1+indgen(i2-i1+1)
nr=n_elements(k)
if nr lt 2 then begin
print,' AVERAGESPEC error: less than 2 files specified - returning'
if keyword_set(stp) then stop,'AVERAGESPEC>>>'
return,0
endif
gdat,file,h,w,f,e,k(0)
time=h(5)
if time lt 0 then time=-60.*time
if time eq 0 then noweight=0
if not keyword_set(noweight) then f=f*time
for j=1,nr-1 do begin
i=k(j)
gdat,file,h,w1,f1,e,i
f1=interpol(f1,w1,w)
if not keyword_set(noweight) then f1=f1*time
f=f+f1
t=h(5)
if t lt 0 then t=-60.*t
time=time+t
endfor
f=f/nr
if not keyword_set(noweight) then f=f/time
if irec then i2=w
print,'Total time in average spectrum = ',time,' seconds'
if keyword_set(save) then begin
if time gt 32767 then time=-fix(time/60.) else time=fix(time)
h(5)=time
h(100)=byte('Average of records '+strtrim(i1,2)+' - '+strtrim(i2,2))
kdat,file,h,w,f,e,-1
endif
if keyword_set(stp) then stop,'AVERAGESPEC>>>'
return,f
end
averagespec.pro;2 000640 000451 000017 00000000000 06615503101 020114 1averagespec.pro ustar 00fwalter users 000002 251444 avspt.pro 000640 000451 000017 00000002506 06615503101 013614 0 ustar 00fwalter users 000002 253260 ;**********************************************************************
function avspt,f,e,kk,bad=bad ;find single point zeros in IUEHI file and remove
np=n_elements(f)-1
if n_elements(bad) gt 0 then k=bad else case 1 of
(min(e) eq 0) and (max(e) lt 10): k=where(e gt 1)
(min(e) lt 0) : k=where(e lt -200)
else : k=where(e eq 0)
endcase
if n_elements(k) le 1 then return,f
nk=n_elements(k)
k1=[k(1:*),np]
k2=[0L,k(0:*)] ;was :np - fixed by FMW 6/24/93
dk1=(k1-k)-1 ; dist to next point
dk2=(k-k2)-1
dk=dk1*dk2
kz=where(dk ne 0,nkz) ;single point zeros
if nkz eq 0 then return,f ;no single point zeros found
kk=k(kz) ;indices of single point zeros
if kk(0) eq 0 then begin
if n_elements(kk) gt 1 then kk=kk(1:*) else return,f ;first point only
endif
nkk=n_elements(kk)
if kk(nkk-1) eq np then begin
if nkk gt 1 then kk=kk(0:nkk-2) else return,f ;last point only
endif
nkk=n_elements(kk)
ff=f
for i=0,nkk-1 do ff(kk(i))=(ff(kk(i)-1)+ff(kk(i)+1))/2.
if (k(0) eq 0) and (k(1) ne 1) then begin
ff(0)=ff(1)
kk=[0,kk]
nkk=nkk+1
endif
if (k(nk-1) eq np) and (k(nk-2) ne np-1) then begin
ff(np)=ff(np-1)
kk=[kk,np]
nkk=nkk+1
endif
if !quiet ne 3 then print,nkk,' single point zeros restored'
return,ff
end
avspt.pro;2 000640 000451 000017 00000000000 06615503101 015632 1avspt.pro ustar 00fwalter users 000002 253260 bargraph.pro 000640 000451 000017 00000001434 06615503101 014251 0 ustar 00fwalter users 000002 251447 ;****************************************
pro bargraph,x,y,ynz=ynz,stp=stp
z=y/max(abs(y))
if keyword_set(ynz) then begin
ymin=min(z) & ymax=max(z)
endif else begin
ymin=-1. & ymax=1.
if min(z) ge 0. then ymin=0.
if max(z) le 0. then ymax=0.
endelse
sp=(ymax-ymin)/20.
nl=22
graph=bytarr(79,nl)+32b
np=n_elements(y)<75
graph(5,*)=byte('|') & graph(78,*)=byte('|')
lab=strarr(nl-1)
for i=0,nl-2 do lab(i)=string(ymin+sp*i,'(F5.2)')
for i=0,nl-2 do graph(0,nl-2-i)=byte(lab(i))
graph(*,0)=byte('_') & graph(*,nl-2)=byte('_')
;
graph(5,nl-1)=byte(strtrim(x(0),2))
graph(5+np-1,nl-1)=byte(strtrim(x(np-1),2))
if sp eq 0. then sp=1.
yp=fix((z-ymin)/sp)
for i=0,np-1 do graph(i+5,nl-2-yp(i))=byte('*')
for i=0,nl-1 do print,string(graph(*,i))
if keyword_set(stp) then stop
return
end
bargraph.pro;2 000640 000451 000017 00000000000 06615503101 016721 1bargraph.pro ustar 00fwalter users 000002 251447 bbody.pro 000640 000451 000017 00000005666 06615503101 013571 0 ustar 00fwalter users 000002 252730 ;***********************************************************************
PRO BBODY,W,F,EPS,temp=temp ; COMPUTE AND OVERPLOT BLACK BODY SPECTRUM
COMMON COM1,H,IK,IFT,NSM,C
COMMON COMXY,Xcur,Ycur,zerr
common icurunits,xunits,yunits,title,col1,col2,col3
IMULT=0 & IPT=0
if n_params(0) lt 2 then return
if n_params(0) lt 3 then eps=fix(w)*0+100
if max(eps) eq 0 then eps=eps+100
E=FLOAT(EPS/(ABS(EPS)>1))
E=E>0.
if n_elements(col2) eq 0 then opcol=!p.color else opcol=col2
REENT:
case 1 of
!d.name eq 'PS' : begin
read,'BBody: scale to all (1), point (2), or range(3), - for multiple: ',itype
if abs(itype) le 1 then begin
x1=min(w) & x2=max(w)
endif
if abs(itype) eq 2 then begin
read,' enter wavelength of point: ',x1
ipt=1
x2=x1
endif
if abs(itype) eq 3 then read,' enter wavelength limits: ',x1,x2
if itype lt 0 then imult=1
I1=FIX(xindex(w,x1)) & I2=FIX(xindex(w,x2))
y1=f(i1)
if ipt eq 1 then i2=0
end
;
else: begin ;plot to screen
print,'BBody: scale to all (0,N), point (1,P), or select range(<space> ,M)'
opstat,' Waiting'
BLOWUP,-1
opstat,' Working'
IF ZERR EQ 90 THEN BEGIN
STOP,' BBODY STOP
GOTO,REENT
ENDIF
IF (zERR EQ 77) OR (ZERR EQ 109) THEN IMULT=1 ;<M>
IF (zERR EQ 80) OR (ZERR EQ 112) THEN BEGIN ;<P>
IMULT=1 & IPT=1
ENDIF
IF zERR EQ 49 THEN IPT=1 ;<1>
IF (zERR EQ 78) OR (ZERR EQ 110) THEN BEGIN ;<N>
IMULT=1 & zERR =48
ENDIF
IF zERR EQ 48 THEN BEGIN ;<0>
I1=0
I2=N_ELEMENTS(W)-1
GOTO,GOON
ENDIF
IF (ZERR EQ 90) OR (ZERR EQ 122) THEN STOP,'BBODY' ;<Z,Z>
IF (ZERR EQ 81) OR (ZERR EQ 113) OR (ZERR EQ 26) THEN RETURN ;<Q,Q,^Z>
X1=Xcur
I1=FIX(xindex(w,x1)) ;TABINV,W,X1,I1
IF IPT EQ 1 THEN BEGIN ;<1,P>
I2=0
Y1=Ycur
ENDIF ELSE BEGIN ;SECOND INPUT
opstat,' Waiting'
BLOWUP,-1
opstat,' Working'
X1=Xcur
I2=FIX(xindex(w,x1)) ; TABINV,W,X1,I2
IF I2 LT I1 THEN BEGIN
T=I2 & I2=I1 & I1=T
ENDIF
IF I2 EQ I1 THEN BEGIN
I2=0
Y1=F(I1)
ENDIF ;ELSE I2=I2-I1+1 ;# POINTS
ENDELSE
end
endcase
GOON:
z='($," Enter BB temperature (K) ")'
if not keyword_set(temp) then begin
print,format=z
READ,TEMP
endif
C1=3.7412E-5 & C2=1.43879
W1=W/1.E8 ;CM
MPLT: BB=C1/W1^5/(EXP(C2/W1/TEMP)-1.)
BB=BB/1.E25 ; NORMALIZE FLUX
FP=F*E
BP=BB*E
IF I2 EQ 0 THEN NORM=Y1/BB(I1) ELSE NORM=TOTAL(FP(I1:I2))/TOTAL(BP(I1:I2))
OPLOT,W,BB*NORM,color=opcol
IK=IK+1
Z='BB:T='+STRTRIM(STRING(TEMP),2)+' N='+STRTRIM(STRING(NORM),2)
print,z
IF IMULT EQ 0 THEN RETURN ELSE BEGIN
z='($,"Enter T, <0 to end ")'
print,format=z
READ,TEMP
IF TEMP LE 0. THEN RETURN
GOTO,MPLT
ENDELSE
END
bbody.pro;2 000640 000451 000017 00000000000 06615503101 015537 1bbody.pro ustar 00fwalter users 000002 252730 bdata.pro 000640 000451 000017 00000001350 06615503102 013527 0 ustar 00fwalter users 000002 252702 ;********************************************************
PRO BDATA,h,X,W,F,E,BW,BF,BDF,V ; SET BAD DATA VECTOR
IF X EQ -1 THEN GOTO,BDAT
WX=X
i=xindex(w,wx)
I=FIX(I+0.5)
IF N_PARAMS(0) LT 8 THEN V=-1111
if (i ge 0) or (i le (n_elements(e)-1)) then E(I)=V
;
BDAT: ; SET BAD DATA VECTORS
if n_elements(e) lt n_elements(f) then e=f*0+100 ;e not passed or scalar passed
if n_elements(h) lt 33 then h33=0 else h33=h(33)
bw=-1 & bf=-1
case 1 of
h33 eq 30: bad=where(e le -1000.,nbad) ; S/N vector
h33 eq 40: bad=where(e lt 0.,nbad) ;error bars
else: bad=where(e lt 0.,nbad)
endcase
;
if nbad gt 0 then begin
BW=W(BAD) & BF=F(BAD)
endif
if n_elements(bdf) eq 1 then bdf=(bdf+1) mod 2
RETURN
END
bdata.pro;2 000640 000451 000017 00000000000 06615503102 015467 1bdata.pro ustar 00fwalter users 000002 252702 bell.pro 000640 000451 000017 00000000421 06615503102 013372 0 ustar 00fwalter users 000002 253262 ;**************************************************************************
pro bell,number,nmax=nmax
if n_elements(nmax) eq 0 then nmax=10
if n_params(0) eq 0 then number=1
for i=1,(number<nmax) do begin
print,FORMAT="($,A1)",string(7b) & wait,0.1
endfor
return
end
bell.pro;2 000640 000451 000017 00000000000 06615503102 015177 1bell.pro ustar 00fwalter users 000002 253262 bigplot.pro 000640 000451 000017 00000003424 06615503102 014124 0 ustar 00fwalter users 000002 252733 ;***************************************************************************
pro bigplot,id,w,f,e,bw,bf ; PLOT MULTIPLE SPECTRAL REGIONS
;vax version
common com1,hd,ik,ift,nsm,c,ndat,ifsm,kblo
common comxy,xcur,ycur,zerr
; read input data
inrec=-9
print,' '
files=''
i1='x'
print,' Enter file, record # for each spectrum (CR,-1 to end)'
while i1 ne '' do begin
read,i1,i2
files=[files,i1]
inrec=[inrec,i2]
endwhile
k=where (files ne '')
files=files(k) & inrec=inrec(k)
k=where (inrec ge 0)
files=files(k) & inrec=inrec(k)
n=n_elements(files)
if n le 0 then return
id=100
READ,' Enter wavelength range',x1,x2
READ,' Enter Ymax', y2
setxy,x1,x2,0.,y2
l=x2-x1
maxbin=l/4096.
; 4096 bins maximum
print,' Minimum bin size=',maxbin,' A'
READ,' Enter bin size (A), or # pts (>1000)',dl
if dl ge 1000 then nbins=fix(dl) else nbins=l/dl
nbins = fix(nbins < 4096)
dl=l/float(nbins)
print,'$I5',nbins,' points to be plotted, dispersion= ', dl
w=x1+dl*findgen(nbins) & f=fltarr(nbins) & e=intarr(nbins)+100
nd=intarr(nbins)
h=intarr(1022)
H(0)=100
h(3)=1000 ;ncam
h(7)=nbins
h(34)=n
hd=h
FUN3,0,0,0,H,HD ;new title
for i=0,n-1 do begin
gdat,files(i),h1,w1,f1,e1,inrec(i)
np=n_elements(w1)
i1=fix(xindex(w,w1(0))) ;first point ;tabinv,w,w1(0),i1
i2=xindex(w,w1(np-1)) ;tabinv,w,w1(np-1),i2
tw=w1(0)+dl*findgen(I2-I1+1) ;section of big vector
rrebin,tw,w1,f1,e1
np=n_elements(f1)
NDT=INTARR(NP)+1
K=WHERE(E1 LT -201,ck)
if ck gt 0 then begin
f1(k)=0. & ndt(k)=0
endif
t=f(i1:i1+np-1)+f1
f(i1)=t
t=nd(i1:i1+np-1)+NDT
nd(i1)=t
t=e(i1:i1+np-1)<e1
e(i1)=t
endfor
f=f/(float(nd)>1.)
K=WHERE((F EQ 0.0) AND (E EQ 100)) ;NO DATA
E(K)=0 ;E=0 IF NO DATA
bdata,hd,-1,w,f,e,bw,bf
pldata,0,w,f,bw,bf,/bdata
return
end
bigplot.pro;2 000640 000451 000017 00000000000 06615503102 016445 1bigplot.pro ustar 00fwalter users 000002 252733 bits.pro 000640 000451 000017 00000000265 06615503102 013424 0 ustar 00fwalter users 000002 253263 ;***********************************************************
pro bits,inp
in=long(inp)
out=intarr(32)
for i=0,30 do out(31-i)=(in and 2L^i)/2L^i
print,'$(4(8I1,1X))',out
return
end
bits.pro;2 000640 000451 000017 00000000000 06615503102 015246 1bits.pro ustar 00fwalter users 000002 253263 blowup.pro 000640 000451 000017 00000006160 06615503103 013766 0 ustar 00fwalter users 000002 260520 ;*****************************************************************
PRO BLOWUP,I,prt=prt,readout=readout ; PROCEDURE TO EXPAND PLOT
; I=-1 FOR CURSOR CALL
; I=0 FOR COMPLETE BLOWUP
; I=1 FOR X ONLY
; I=2 FOR Y ONLY
COMMON COMXY,XCUR,YCUR,ZERR
common vars,var1,var2,var3,var4,var5,psdel,prffit,vrot2
if n_elements(xcur) eq 0 then xcur=mean(!x.crange)
if n_elements(ycur) eq 0 then ycur=mean(!y.crange)
;if n_elements(x) eq 0 then x=200
;if n_elements(y) eq 0 then y=200
;if n_elements(var3) eq 0 then var3=0
;NAX=rdbit(var3,0)
FORM1="($,A7)"
if strupcase(!d.name) eq 'X' or strupcase(!d.name) eq 'WIN' then isrc=1 else isrc=0
IF I EQ 91 THEN GOTO,EXPAND ;EXPAND PLOT ;[
IF I GE 30 THEN GOTO,SHFT ;<,>,^,6
IF I EQ -1 THEN BEGIN
if (xcur le !x.crange(0)) or (xcur gt !x.crange(1)) then xcur=mean(!x.crange)
if (ycur le !y.crange(0)) or (ycur gt !y.crange(1)) then ycur=mean(!y.crange)
case 1 of
isrc: begin
tvcrs,xcur,ycur,/data
if n_elements(readout) ne 0 then kwhere,readout else begin
key=" "
key=get_kbrd(1)
zerr=FIX(byte(key))
ZERR=ZERR(0)
CURSOR,xcur,ycur,/NOWAIT
endelse
end
strupcase(!d.name) eq 'TEK': begin
CURSOR,xcur,ycur
zerr=!err
key=" "
key=get_kbrd(1)
end
ELSE: begin
key=get_kbrd(1)
zerr=FIX(byte(key))
ZERR=ZERR(0)
end
endcase
if keyword_set(prt) then print,xcur,ycur
RETURN
ENDIF ;i=-1
;
; second call for i=0,1,2
;
XMAX=!X.range(1)
XMIN=!X.range(0)
YMAX=!Y.range(1)
YMIN=!Y.range(0)
X1=Xcur
Y1=Ycur
case 1 of
isrc: begin
OPSTAT,' Waiting'
tvcrs,xcur,ycur,/data
z0=zerr
if n_elements(readout) ne 0 then kwhere,readout else begin
key=" "
key=get_kbrd(1)
CURSOR,xcur,ycur,/NOWAIT
endelse
zerr=z0
end
strupcase(!d.name) eq 'TEK': begin
PRINT,FORMAT=FORM1,'Waiting'
CURSOR,xcur,ycur,/WAIT,/data
zerr=!err
key=" "
key=get_kbrd(1)
PRINT,FORMAT=FORM1,' '
print,format="($,A)",string("12b)
end
ELSE: begin
key=" "
key=get_kbrd(1)
end
endcase
X2=Xcur
Y2=Ycur
!x.range=[(x1<x2),(x1>x2)]
!y.range=[(y1<y2),(y1>y2)]
;
IF I EQ 1 THEN !Y.range=[YMIN,ymax] ; EXPAND X AXIS ONLY
IF I EQ 2 THEN !X.range=[XMIN,xmax] ; EXPAND Y AXIS ONLY
ZERR=68
RETURN
;
EXPAND: ;EXPAND SCALES 2 TIMES
RX=(!X.range(1)-!X.range(0))
CX=RX/2.+!X.range(0)
RY=(!Y.range(1)-!Y.range(0))
CY=RY/2.+!Y.range(0)
!X.RANGE=[CX-RX,CX+RX]
!Y.RANGE=[CY-RY,CY+RY]
I=68
RETURN
;
SHFT: ; SHIFT ALONG AXES
; > (62) SHIFT UP 1 FRAME (X)
; < (60) SHIFT DOWN 1 FRAME (X)
; ^ (94) SHIFT UP IN Y
; 6 (54) SHIFT DOWN IN Y
DX=!X.crange(1)-!X.crange(0)
DY=!Y.crange(1)-!Y.crange(0)
ISIGN=1
IF (I EQ 60) OR (I EQ 62) THEN BEGIN
IF I EQ 60 THEN ISIGN=-1
!X.range(0)=!X.crange(0)+ISIGN*DX
!X.range(1)=!X.crange(1)+ISIGN*DX
ENDIF
IF (I EQ 94) OR (I EQ 54) THEN BEGIN
IF I EQ 54 THEN ISIGN=-1
!Y.range(0)=!Y.crange(0)+ISIGN*DY
!Y.range(1)=!Y.crange(1)+ISIGN*DY
ENDIF
I=68
RETURN
END
blowup.pro;2 000640 000451 000017 00000000000 06615503103 016157 1blowup.pro ustar 00fwalter users 000002 260520 centrd.pro 000640 000451 000017 00000001514 06615503103 013744 0 ustar 00fwalter users 000002 252735 ;*************************************************************
PRO CENTRD,WAVE,FLUX,YD,I1,I2,XCEN,FD ;CENTROID
; input I1,I2 in bins
; output xcen in units of wave
XCEN=-1.
IF I2 LT I1 THEN BEGIN
T=I2
I2=I1
I1=T
ENDIF
i1=(i1>0) ;force I1 > 0
np=n_elements(wave)-1
i2=i2<np ;force i2 < greatest element
IF (I2-I1) LE 1. THEN RETURN
FD=FLUX-YD
IF TOTAL(FD(i1+1:(i2-1)>(I1+1))) LT 0 THEN FD=-FD
Q=WAVE
P=WAVE*FD
; THE FOLLOWING CODE IS MODIFIED FROM [177,1]INTEG
ILO=LONG(I1+1)
IHI=LONG(I2)
; COMPUTE SUM OF WAVE*FD
xcen=total(p(ilo:ihi))
xcen=xcen+p(ilo-1)*(float(ilo)-i1)+p((ihi+1)<np)*(i2-float(ihi))
; NOW COMPUTE SUM OF FLUX VECTOR
P=FD
xc=total(p(ilo:ihi))
xc=xc+p(ilo-1)*(float(ilo)-i1)+p((ihi+1)<np)*(i2-float(ihi))
XCEN=XCEN/XC ; CENTROID IS SUM WAVE*FD / SUM OF FD
fd=flux(i1:i2)
RETURN
END
centrd.pro;2 000640 000451 000017 00000000000 06615503103 016106 1centrd.pro ustar 00fwalter users 000002 252735 check_stdfile.pro 000640 000451 000017 00000001074 06615503104 015251 0 ustar 00fwalter users 000002 260524 ;**************************************************************
pro check_stdfile
common icdisk,icurdisk,icurdata,ismdata,objfile,stdfile,userdata,recno
;
if stdfile eq 'nofile' then begin
read,' standard file is not defined. Please enter name here: ',stdfile
endif
;
if not ffile(stdfile+'.icd') then case 1 of
ffile(userdata+stdfile+'.icd'): stdfile=userdata+stdfile
ffile(icurdata+stdfile+'.icd'): stdfile=icurdata+stdfile
else: begin
print,'File ',stdfile,' not found - returning
stdfile='nofile'
; return
end
endcase
return
end
check_stdfile.pro;2 000640 000451 000017 00000000000 06615503104 020722 1check_stdfile.pro ustar 00fwalter users 000002 260524 checklim.pro 000640 000451 000017 00000001667 06615503104 014257 0 ustar 00fwalter users 000002 252736 ;*****************************************************************
pro checklim,f,lam1,lam2,e
common comxy,xcur,ycur,zerr
if n_params(0) lt 3 then return
if !d.name eq 'PS' then return
np=n_elements(f)
w=findgen(np)
lam0=[lam1,lam2]
psym=!p.psym
!x.title='bin'
!y.title='flux'
!p.title='CHECKLIM'
setxy
!p.position=[.2,.2,.9,.9]
;
restart:
!c=-1
plot,w,f,psym=0
xcur=mean(!x.crange) & ycur=mean(y.crange)
if n_params(0) ge 4 then begin
k=wherebad(e,1)
oplot,w(k),f(k)*0.+!y.crange(0),psym=1
endif
oplot,[lam1,lam1],!y.crange
oplot,[np-lam2+1,np-lam2+1],!y.crange
z=' Use cursor to mark wavelength region, 0 if OK, 2 to restart'
print,z
blowup,-1
if zerr eq 48 then goto,ret ;return
if zerr eq 50 then goto,restart
lam1=long(xcur)>20L
blowup,-1
lam2=long(xcur)
if lam2 lt lam1 then begin
t=lam2
lam2=lam1
lam1=t
endif
lam2=(np+1-lam2)>20L
print,lam1,lam2
setxy,lam1,np+1-lam2
goto,restart
;
ret:
setxy
return
end
checklim.pro;2 000640 000451 000017 00000000000 06615503104 016710 1checklim.pro ustar 00fwalter users 000002 252736 chisq.pro 000640 000451 000017 00000000552 06615503104 013572 0 ustar 00fwalter users 000002 260525 ;*************************************************************************
function chisq,y,dy,model
; return chi-square value
if n_params(0) eq 0 then return,-1
if n_params(0) eq 1 then dy=sqrt(y)
if n_params(0) lt 3 then model=mean(y)
n=n_elements(y)
ddy=dy & k=where(ddy eq 0.,nk) & if nk gt 0 then ddy(k)=1.
z=(y-model)/ddy
chisq=total(z*z)
return,chisq
end
chisq.pro;2 000640 000451 000017 00000000000 06615503104 015563 1chisq.pro ustar 00fwalter users 000002 260525 chlo.pro 000640 000451 000017 00000003477 06615503104 013422 0 ustar 00fwalter users 000002 260526 ;***************************************************
pro chlo,infile,outfile,r1,r2 ;change iue low dispersion format 0 to format 8
if n_params(0) eq 1 then begin
if ifstring(infile) ne 1 then begin ;not a string
if infile eq -1 then begin
print,' CHLO : change ICUR data format 0 to 8'
print,' call CHLO,infile,outfile,r1,r2'
print,' infile: input file name (format 0), default=LODAT'
print,' outfile: output file name (format 8), default=IUELO'
print,' r1: first record to transcribe, default=0'
print,' r2: last record, =r1 if r1 specified, 999 otherwise'
return
endif
endif
endif
if n_params(0) eq 0 then infile='LODAT'
infile=infile+'.DAT'
if n_params(0) lt 2 then outfile='IUELO'
if n_params(0) lt 3 then begin
r1=0 & r2=999
endif
if n_params(0) eq 3 then r2=r1 ;only one record
r1=(r1>0)
r2=(r2>r1)
on_ioerror,nofile
get_lun,lu
openr,lu,outfile+'.dat'
close,lu
goto,proceed
;
nofile: openw,lu,outfile+'.dat/recs:1000',8196
;
proceed: close,lu
free_lun,lu
;
ws=1100.+findgen(900) ;SWP wavelength vector
wl=1900.+1.5*findgen(934) ;LWP/R wavelength vector
FOR I=r1,r2 do BEGIN
on_ioerror,done
GDAT,infile,h,w,f,e,i
if n_elements(h) eq 1 then goto,done
if n_elements(h) gt 1020 then h=h(0:1019)
ncam=h(3)
print,'record=',i,' NCAM=',ncam
if ncam le 2 then begin ;LWP/LWR
w0=wl
dl=1.5
endif else begin
w0=ws
dl=1.0
endelse
f1=interpol(f,w,w0)
e1=fix(interpol(e,w,w0))
k=where(w0 le max(w)) ;valid points
w0=w0(k)
f1=f1(k)
e1=e1(k)
h(20)=fix(w0(0))
h(21)=1000.*(w0(0)-float(h(20)))
h(22)=fix(dl)
h(23)=1000.*(dl-float(h(22)))
; h(100:159)=32b
h(0)=8
kdat,outfile,h,w0,f1,e1,i
endfor
done:
return
end
chlo.pro;2 000640 000451 000017 00000000000 06615503104 015220 1chlo.pro ustar 00fwalter users 000002 260526 chtitle.pro 000640 000451 000017 00000001432 06615503105 014112 0 ustar 00fwalter users 000002 260530 ;************************************************************************
pro chtitle,file,rec,title,helpme=helpme
if n_params(0) eq 0 then helpme=1
if ifstring(file) ne 1 then helpme=1 ;integer passed
if keyword_set(helpme) then begin ;help
print,' '
print,'* CHTITLE - change title in header vector'
print,'* calling sequence: CHTITLE,file,rec,[title]'
print,'* file: name of .ICD data file'
print,'* rec: record number to be updated'
print,'* title: optional string title (prompted if not present)'
print,'* both arguments are prompted for if not passed'
print,' '
return
endif
;
if n_params(0) lt 2 then read,' Enter record number: ',rec
gdat,file,h,w,f,e,rec
ctit,h,nch,title
if nch ne 1 then kdat,file,h,w,f,e,rec
return
end
chtitle.pro;2 000640 000451 000017 00000000000 06615503105 016432 1chtitle.pro ustar 00fwalter users 000002 260530 cleanspec.pro 000640 000451 000017 00000005077 06615503105 014425 0 ustar 00fwalter users 000002 260531 ;***********************************************
pro cleanspec,file,rec,helpme=helpme,stp=stp,autosave=autosave
common comxy,xcur,ycur,zerr
common com1,H,IK,IFT,NSM,C,NDAT,ifsm,kblo,h2,ipdv,ihcdev
if keyword_set(helpme) then begin
print,' '
print,'* CLEANSPEC - clean bad pixels in .ICD files'
print,'* calling sequence: CLEANSPEC,file,rec'
print,'* FILE: name of .ICD file'
print,'* REC: record to clean, default=all'
print,'* '
print,'* KEYWORDS:'
print,'* AUTOSAVE: automatically save cleaned spectrum; inquire if not set'
print,'* '
print,' '
return
endif
;
if n_elements(rec) eq 0 then rec=-999
if rec(0) eq -999 then begin
icdfile,file,rec
rec=indgen(rec+1)
endif
nrec=n_elements(rec)
p0=!p.psym
gc,11
setxy
for irec=0,nrec-1 do begin
gdat,file,h,w,f,e,rec(irec)
if (n_elements(h) lt 2) and (h(0) eq -1) then return
!p.title=get_title(h)
etype=h(33)
plot,w,f,/ynoz
if !d.name eq 'X' then wshow
zerr=1
change=0
while zerr eq 1 do begin ;do a spectrum
print,' click left button on bad point, right to quit'
cursor,xcur,ycur
zerr=!err
if zerr eq 1 then begin
change=1
!p.psym=10
locate,-1,w,f,nbins=25
print,' mark left and right points to interpolate between'
cursor,xl,y,wait=3
print,' mark other point'
wait,0.2
cursor,xr,y,wait=3
i1=fix(xindex(w,xl)+0.5)
i2=fix(xindex(w,xr)+0.5)
IF I2 LT I1 THEN BEGIN
T=I1 & I1=I2 & I2=T
T=Xl & Xl=Xr & xr=t
ENDIF
nb=i2-i1
print,i1,i2
if nb gt 0 then begin
SL=(F(I2)-F(I1))/FLOAT(NB)
FOR I=I1+1,I2-1 DO F(I)=F(I1)+SL*FLOAT(I-I1)
case 1 of
etype eq 0:
etype eq 1:
etype eq 10: FOR I=I1+1,I2-1 DO e(I)=0
etype eq 20: FOR I=I1+1,I2-1 DO e(I)=-200
etype eq 30: FOR I=I1+1,I2-1 DO e(I)=-e(I1)+SL*FLOAT(I-I1)
else:
endcase
endif
oplot,w,f,color=15
wait,1. & setxy & plot,w,f,/ynoz
endif else zerr=0 ; interpolate one line
endwhile ; do a spectrum
; setxy
; plot,w,f
if (change) then begin
if keyword_set(autosave) then igo='0' else begin
igo=''
read,' OK? enter 0 to save, anything else to go on: ',igo
endelse
if strtrim(igo,2) eq '0' then kdat,file,h,w,f,e,rec(irec)
endif
endfor
!p.psym=p0
if keyword_set(stp) then stop,'CLEANSPEC>>>'
return
end
cleanspec.pro;2 000640 000451 000017 00000000000 06615503105 017235 1cleanspec.pro ustar 00fwalter users 000002 260531 closeall.pro 000640 000451 000017 00000000207 06615503105 014257 0 ustar 00fwalter users 000002 260534 ;*****************************************************************************
pro closeall,dum
for i=100,128 do free_lun,i
return
end
closeall.pro;2 000640 000451 000017 00000000000 06615503105 016742 1closeall.pro ustar 00fwalter users 000002 260534 coadd.pro 000640 000451 000017 00000005362 06615503105 013543 0 ustar 00fwalter users 000002 260535 ;**********************************************************************
PRO COADD,W,F,E,nbin,error=error,wave=wave,helpme=helpme
;COADD DATA INTO COARSER BINS
common com1,h,IK,IFT,NSM,C,NDAT,IFSM,KBLO,H1,ipdv,ihcdev
npar=n_params(0)
if npar eq 0 then helpme=1
if npar eq 1 then begin
if n_elements(w) le 2 then helpme=1
endif
if keyword_set(helpme) then begin
print,' '
print,'* COADD - coadd bins together'
print,'* calling sequence: COADD,W,F,E,NBINS, or COADD,W,F,NBINS,'
print,'* or COADD,W,NBINS, or COADD,W'
print,'* W,F,E : 1-3 variables. if 2 or 3 are passed, first is treated'
print,'* as a wavelength or a time'
print,'* a mean F is returned; E returns minimum of E'
print,'* NBINS: if the last parameter is an integer, it is the'
print,'* number of bins to coadd, otherwise input will be requested'
print,'*'
print,'* KEYWORDS:'
print,'* ERROR: if set, W is treated as errors and added in quadrature'
print,'* WAVE: if set, W is treated as a wavelength vector'
print,' '
return
endif
;
if n_elements(h) lt 33 then h33=0 else h33=h(33)
if n_params(0) lt 4 then nbin=-1 else nbin=abs(fix(nbin))>1
case 1 of
npar eq 2: if n_elements(f) eq 1 then nbin=abs(fix(f))>1
npar eq 3: if n_elements(e) eq 1 then nbin=abs(fix(e))>1
else:
endcase
if nbin ne -1 then nvars=npar-1 else nvars=npar ;number of variables passed
;
if nbin eq -1 then begin
PRINT,' '
READ,' HOW MANY BINS TO COADD? ',NBIN
endif
IF NBIN LE 1 THEN RETURN
;
np=n_elements(w)
if np le 1 then np=n_elements(f)
S=np/NBIN
tw=w
;
case 1 of
keyword_set(error): begin ;vector is SNR vector
w=fltarr(s)
ii=nbin*indgen(s)
for i=0,nbin-1 do w=w+tw(ii+i)*tw(ii+i)
w=sqrt(w)/float(nbin)
end
;
keyword_set(wave): begin ;wavelength or time vector
i=indgen(s)*nbin
w=tw(i)
end
;
nvars eq 1: begin ;only one vector passed
w=fltarr(s)
ii=nbin*indgen(s)
for i=0,nbin-1 do w=w+tw(ii+i)
w=w/float(nbin)
end
;
(nvars gt 1) and (n_elements(w) gt 1): begin ;wavelength or time vector
i=indgen(s)*nbin
w=tw(i)
end
nvars ge 2: begin
tf=f
F=FLTARR(S)
ii=nbin*indgen(s)
for i=0,nbin-1 do f=f+tf(ii+i)
f=f/float(nbin)
end
(nvars eq 3) and (n_elements(e) gt 1): begin
if (h33 eq 30) or (h33 eq 40) then te=e*e else te=e
e=FLTARR(S)
ii=nbin*indgen(s)
for i=0,nbin-1 do e=e+te(ii+i)
case 1 of
h33 eq 30: e=sqrt(e) ;add S/N vectors
h33 eq 40: e=sqrt(e)/float(nbin)
else:
endcase
end
endcase
;
if n_elements(h) gt 52 then h(53)=nbin
RETURN
END
coadd.pro;2 000640 000451 000017 00000000000 06615503105 015473 1coadd.pro ustar 00fwalter users 000002 260535 combine.pro 000640 000451 000017 00000011077 06615503105 014101 0 ustar 00fwalter users 000002 260540 ;****************************************************************************
pro combine,infile,inrec,outfile,outrec,average=average,noscale=noscale, $
sumup=sumup,examine=examine,helpme=helpme,stp=stp,edit=edit, $
nosave=nosave ;,etype=etype,sigma=sigma
if keyword_set(helpme) then begin
print,' '
print,'* COMBINE - combine 3 or more spectra with median filtering'
print,'* averages 2 spectra'
print,'* calling sequence: COMBINE,infile,inrec,outfile,outrec'
print,'* INFILE: input .ICD file'
print,'* INREC: array of record numbers to combine'
print,'* OUTFILE: output .ICD file, default=INFILE'
print,'* OUTREC: output record, default=next free record'
print,'*'
print,'* KEYWORDS:'
print,'* AVERAGE: use averaging, not median filtering'
; print,'* ETYPE: set to force ETYPE when .ICD file is created'
print,'* EXAMINE: examine data at all stages.'
print,'* NOSCALE: do not scale flux vectors'
; print,'* SIGMA: sigma limit for median filtering
print,'* SUMUP: sum all flux, no scaling - weights by exp. time'
print,' '
return
endif
;
if keyword_set(edit) then examine=1
if n_elements(infile) eq 0 then begin
infile=''
read,' Enter name of .ICD file: ',infile
endif
if n_elements(inrec) eq 0 then begin
inrec=-1 & ii=0
print,' Enter list of record numbers, -1 to end'
while ii ge 0 do begin
read,'>>',ii
inrec=[inrec,ii]
endwhile
k=where(inrec ge 0,nk)
if nk lt 2 then begin
print,' Fewer than 2 records specified - returning'
return
endif
inrec=inrec(k)
endif
if n_elements(outfile) eq 0 then outfile=infile
if n_elements(outrec) eq 0 then outrec=-1
nrec=n_elements(inrec)
if nrec lt 2 then begin
print,' Fewer than 2 records specified - returning'
return
endif
;
gdat,infile,h,w,f,e,inrec(0)
if keyword_set(examine) then plot,w,f,/ynoz
ft=f*h(5) ;flux
np=h(7)
kp=indgen(np/2)+np/4 ;middle half of vector for scaling
tf=total(f(kp))
h0=h
h0(9)=1
epstype=h(33)
darr=fltarr(np,nrec) & earr=darr
darr(0,0)=f & earr(0,0)=e
scalef=fltarr(nrec)+1.
;
for i=1,nrec-1 do begin
gdat,infile,h,w,f,e,inrec(i)
ft=ft+f*h(5) ;total flux
scalef(i)=tf/total(f(kp))
darr(0,i)=f & earr(0,i)=e
h0(5)=h0(5)+h(5) ;increment time
h0(9)=h0(9)+1
if keyword_set(examine) then begin
print,' Record',i,': scale factor= ',scalef(i),' Accum time= ',h0(5)
if not keyword_set(noscale) then sf=scalef(i) else sf=1.
oplot,w,f*sf
endif
endfor
if keyword_set(edit) then stop,' Edit data now'
;
if (epstype eq 0) and (not keyword_set(sumup)) then average=1
if (nrec eq 2) and (not keyword_set(sumup)) then average=1
;
if keyword_set(noscale) then scale=0 else scale=1 ;scaling turned off
;
if scale eq 1 then for i=1,nrec-1 do begin ;scale flux vector
f=darr(*,i)*scalef(i)
darr(0,i)=f
endfor
fbar=sum(darr,1)/nrec ;average flux
;
if keyword_set(examine) then begin
if not keyword_set(sumup) then sumup=0
if not keyword_set(average) then average=0
if not keyword_set(noscale) then noscale=0
print,'Average:',average,' Sumup:',sumup,' Noscale:',noscale,' Etype:',epstype
endif
;
case 1 of
keyword_set(sumup): f=ft/h0(5) ;total flux
keyword_set(average): f=fbar ;average flux
else: begin ;median filter
f=medarr(darr) ;median array
; if not keyword_set(sigma) then sigma=5.
; csig3,sigma,nrec,darr,earr,f,expo
;k=where(expo eq 0,nzero)
; carr=narr/(carr > .001) ;sigma
; narr=narr/float((expo > 1)) ;counts per second
; if nzero gt 0 then begin ;no exposure - average adjacent points
; print,'nzero,element =',nzero,k
; for i=0,nzero-1 do begin
; narr(k(i))=(narr(k(i)-1)+narr(k(i)+1))/2.
; carr(k(i))=(carr(k(i)-1)+carr(k(i)+1))/2.
end
endcase
;
case 1 of ;coadd errors
epstype eq 30: begin
e=sqrt(sum(earr*earr,1)) ;s/n vector
end
else: begin
e=fltarr(np)
for i=0,np-1 do e(i)=f(i)/(stddev(darr(i,*))>1.e-25) ;standard deviation
h0(33)=30
end
endcase
;
if keyword_set(examine) then begin
gc,11
oplot,w,f,color=55
erbar,2,f,f/e,w,color=15
wshow
; stp=1
endif
;
if not keyword_set(nosave) then kdat,outfile,h0,w,f,e,outrec ;,epstype=etype
;
if keyword_set(stp) then stop,'COMBINE>>>'
return
end
combine.pro;2 000640 000451 000017 00000000000 06615503105 016373 1combine.pro ustar 00fwalter users 000002 260540 contfl.pro 000640 000451 000017 00000011242 06615503106 013746 0 ustar 00fwalter users 000002 260550 ;**************************************************************************
PRO CONTFL,file,RECS,wmin,wmax,fluxes,out=out,helpme=helpme,stp=stp,plt=plt, $
hcpy=hcpy,abdor=abdor,title=title,two=two
;
if n_params(0) lt 2 then helpme=1
if keyword_set(helpme) then begin
print,' '
print,'* CONTFL - measure continuum fluxes'
print,'* calling sequence: CONTFL,FILE,RECS,wmin,wmax,f'
print,'* FILE: name of data file'
print,'* RECS: record numbers to be measureD, -1 for all'
print,'* WMIN,WMAX: bandpass(es) to be integrated'
print,'* F: output integrated fluxes'
print,'*'
print,'* KEYWORDS:'
print,'* OUT: set to print listing'
print,' '
return
endif
;
if not ifstring(file) then begin
print,' file must be a string'
return
endif
S=n_elements(RECS)
IF S EQ 1 THEN BEGIN ;1 ENTRY
T=RECS
RECS=INTARR(1)+T
ENDIF
if recs(0) Le -1 then recs=INDGEN(GET_NSPEC(FILE))
N=N_ELEMENTS(RECS)
PRINT,N,' RECORDS TO BE MEASURED'
fluxes=fltarr(n,3)
jd=dblarr(n)
;
if n_params(0) eq 2 then begin
WMIN=-1.
WMAX=-1.
READ,' ENTER PAIRS OF WAVELENGTHS, -1 TO END'
READ,WMIN,WMAX
IF WMIN EQ -2 THEN BEGIN ;swp continuua
WMIN=[1350.,1450.,1570.,1680.]
WMAX=[1390.,1500.,1620.,1730.]
GOTO,BAIL
ENDIF
if WMIN eq -3 then begin ;MgII basal fluxes
wmin=[2789.5,2793.5,2800.5,2803.]
wmax=wmin+4.
endif
IF WMIN LT 0. THEN GOTO, BAIL
while wmin gt 0. do begin
read,w1,w2
wmin=[wmin,w1]
wmax=[wmax,w2]
endwhile
endif
BAIL:
if n_elements(wmin) eq 1 then begin
wmin=wmin+fltarr(1)
wmax=wmax+fltarr(1)
endif
;
ffmt='(F9.3)'
tfmt='(E11.3)'
if keyword_set(out) then BEGIN
IF IFSTRING(out) then ofile=out else ofile='contfl'
if strlen(get_ext(ofile)) eq 0 then ofile=ofile+'.lst'
OPENW,lu,ofile,/get_lun
endif else lu=-1
PRINTF,lu,' CONTINUUM FLUX MEASUREMENTS'
printf,lu,' Data from ',file
L=WHERE (WMIN GT 0.)
L=MAX(L)
PRINTF,lu,' The ',L+1,' wavelengths ranges are:'
for i=0,l do printf,lu,string(wmin(i),ffmt),string(wmax(i),ffmt)
printf,lu,' '
FOR I=0,N-1 DO BEGIN
GDAT,file,H,W,F,E,recs(i)
IF N_ELEMENTS(H) LT 2 THEN GOTO,DONE
juldate,[1900+h(12),h(10),h(11),h(13),h(14)],jdy
jd(i)=jdy
IF h(3) le 4 THEN BEGIN
CAMERA=STRMID(' LWP LWR SWP SWR',H(3)*4,4)
image=h(4)
if image lt 0 then image=image+65636L
PRINTF,lu,'Record',RECS(I),' IUE camera=',H(3),' image=',image
endif else BEGIN
printf,lu,'Record',RECS(I),' NCAM=',string(H(3),'(I4)'), $
': ',STRtrim(BYTE(H(100:159)>32),2)
endelse
FOR J=0,L DO BEGIN
IF ((WMIN(J) GT MAX(W)) OR (WMAX(J) LT W(0))) THEN GOTO,SKIP
i1=xindex(w,wmin(j)) ;TABINV,W,WMIN(J),I1
i2=xindex(W,WMAX(J))
I1=FIX(I1+0.5) & I2=FIX(I2+0.5)
I2=I2+1
TF=TOTAL(F(I1:I2))*(WMAX(J)-WMIN(J))/(I2-I1)
k=i1+indgen(i2+1-I1)
TFA=TF/(WMAX(J)-WMIN(J))
BADE=E(I1:I2)
NBADE=WHERE(BADE LE -1600,count)
IF count gt 0 THEN NB=N_ELEMENTS(NBADE) ELSE NB=0
z=string(WMIN(J),ffmt)+string(WMAX(J),ffmt)
Z=Z+' Flux,Flux/A='+string(TF,tfmt)+' '+string(TFA,tfmt)+' '
case h(33) of
30: begin
k0=where(e eq 0,nk)
if nk gt 0 then begin
e(k0)=1.
eb=f/e
for ik=0,nk-1 do eb(k0)=(eb((k0-1)>0)+eb(k0+1))/2.
endif else eb=f/e
tfe=total(sqrt(eb(k)*eb(k)))/(n_elements(k)-1)
end
40: tfe=total(sqrt(e(k)*e(k)))/(n_elements(k)-1)
else: tfe = 0.
endcase
if j eq 0 then begin
fluxes(i,0)=tf & fluxes(i,1)=tfa & fluxes(i,2)=tfe
endif
if h(3) le 4 then $
z=z+STRING(NB,'(I3)')+' points of'+STRING(I2,'(I4)')+' are saturated'
PRINTF,lu,'Range:',z
SKIP:
ENDFOR
ENDFOR
DONE:
;
if n_elements(hcpy) gt 0 then plt=1
if keyword_set(plt) then begin
x=indgen(n)
!x.title='!6 index'
ztitle='!6 CONTFL'
y=fluxes(*,0)
eb=fluxes(*,2)
if keyword_set(abdor) then begin
phase0=2444296.575D0 & period=0.51479D0
jd=jd+2400000.0D0
cycle=(jd-phase0)/period
x=cycle-long(cycle(0))
!x.title='!6 phase'
ztitle='!6 AB Dor'
if keyword_set(two) then begin
x=[x,x+1.] & y=[y,y] & eb=[eb,eb]
endif
endif
if n_elements(title) eq 0 then !p.title=ztitle else !p.title=title
if n_elements(hcpy) gt 0 then sp,'ps'
plot,x,y,/ynozero
erbar,2,y,eb,x
if !d.name eq 'X' then wshow
make_hcpy,hcpy
endif
;
if keyword_set(out) then begin
CLOSE,lu
free_lun,lu
print,' Listing is in ',ofile
endif
if keyword_set(stp) then stop,'CONTFL>>>'
RETURN
END
contfl.pro;2 000640 000451 000017 00000000000 06615503106 016117 1contfl.pro ustar 00fwalter users 000002 260550 copy_icur.tex 000640 000451 000017 00000005346 06615503106 014471 0 ustar 00fwalter users 000002 260545 %*******************************************************************
\documentstyle[12pt]{article}
\newcommand{\skp}{\mbox{ }\\}
\oddsidemargin=0pt
\headheight -40pt
\textheight=9.00in
\textwidth=6.75in
\pagestyle{empty}
\begin{document}
\begin{center}{\bf How to copy ICUR to your home institution}\end{center}
All files are in {\it 44156::rulupi\$dka200:[users.public.icur]}. Note that
this is a VMS machine. There are two ways to access the files.
\begin{enumerate}
\item not-so-anonymous FTP.
\begin{verbatim}
ftp sbast1.ess.sunysb.edu
>logon as guest, password guest
cd rulupi$dka200:[users.public.icur]
ls (if you want to know what's there)
get ...
(Use binary mode for the .TBL and .ICD files; all other files are ASCII.)
\end{verbatim}
\item VMS copy
\begin{verbatim}
dir 44156"guest guest"::rulupi$dka200:[users.public.icur]*.*
(to see what's there)
copy 44156"guest guest"::rulupi$dka200:[users.public.icur]xxx.xxx *
\end{verbatim}
Note that the guest account is restricted - you cannot log into it to do
anything. You can do remote copies and directories, however.
\end{enumerate}
The files you need are:
\begin{description}
\item[ICUR.COM] a sample VMS command file for running ICUR.
\item[ICUR.TLB] the VMS text library containing all the procedures
(1461 blocks). UNIX users can grab
all the individual procedures from the [.pro] subdirectory.
\item[ICUR.TEX] the LATEX documentation file (42 pages).
\item[ICURSTARTUP\_V2.PRO] a main procedure which runs GOICUR.
\item[ICUR.MSG] an optional file typed by GOICUR, containing messages.
\item[*.LIN] 2 files (UV and OPT) containing incomplete line lists.
\item[*.DAT] the 5 files
NANDY.DAT, ORI.DAT, SAVMAT.DAT, SEATON.DAT, and SMC.DAT,
contain interstellar extinction curves. These are
needed for the u command. There are also
12 PSF*.DAT files, useful only for persons having pre-COSTAR
GHRS
spectra and access to the GHRS software, for doing
quick-and-dirty deconvolutions.
\item[SAMPLE.ICD] sample data file.
\item[EXPORT.TXT] the list of current updates to the package (a feeble attempt
at configuration control).
\end{description}
Note for UNIX users: The files may still be available via anonymous FTP from
ftp.astro.psu.edu, in directory pub/nefftp/icur The file
README contains all installation instructions.
The distribution is contained in the compressed TAR file icur28may92.tar.z
This tar file is currently out of date, and its use is not recommended.
Please report any problems to fwalter@astro.sunysb.edu.
\end{document}
copy_icur.tex;2 000640 000451 000017 00000000000 06615503106 017341 1copy_icur.tex ustar 00fwalter users 000002 260545 corrflux.pro 000640 000451 000017 00000007427 06615503106 014340 0 ustar 00fwalter users 000002 260560 ;***********************************************************************
pro corrflux,file,rec,w0def=w0def,nosave=nosave,ffx=ffx,helpme=helpme, $
stp=stp,dw=dw,noplot=noplot,ffl=ffl,nsm=nsm,ddw=ddw
common comxy,xcur,ycur,zerr
if not keyword_set(w0def) then w0def=6563.
if n_params(0) eq 0 then helpme=1
if n_elements(file) eq 0 then file='-1'
if strtrim(string(file),2) eq '-1' then helpme=1
hlp:
if keyword_set(helpme) then begin
print,' '
print,'* CORRFLUX - procedure to tweak up flux corrections, using correction'
print,'* generated by FUDGEFLUX.'
print,'*'
print,'* Calling sequence: CORRFLUX,file,recs'
print,'* file: name of data file, 6 = ECHEL.DAT, -1 for this help message.'
print,'* recs: -1 for this help message'
print,'* an integer or an array of record numbers to be done interactively'
print,'* ''*'' for automatic correction of all'
print,'* KEYWORDS:
print,'* ffx: name of .FFX file, 0 for FUDGE, containing correction factor'
print,'* ffl: if set, use .FFL file rather than .FFX file
print,'* w0def: the wavelength of the plot region; default=',strtrim(w0def,2)
print,'* NOSAVE: set to see results without overwriting data file.'
print,' '
return
endif
;
case 1 of
n_elements(ffl) gt 0: begin
ffx=file
if not ffile(ffx+'.ffl') then ffx=strtrim(ffl,2)
if not ffile(ffx+'.ffl') then ffx='fudge'
if not ffile(ffx+'.ffl') then begin
print,' correction factor file ',ffx+'.FFL not found. Returning'
return
endif else wl=get_ffl(ffx) ;correction factor
nl=n_elements(wl)
fl=wl*0.
ffl=1
end
else: begin
if not keyword_set(ffx) then ffx=file
if not ffile(ffx+'.ffx') then ffx='fudge'
if not ffile(ffx+'.ffx') then begin
print,' correction factor file ',ffx+'.FFX not found. Returning'
return
endif else fact0=get_ffx(ffx,w0) ;correction factor
ffl=0
end
endcase
;
if n_elements(rec) eq 0 then $
read,' enter record number, -1 for help,''*'' for all: ',rec
if ifstring(rec) eq 1 then begin ;string passed
if rec ne '*' then begin
print,' CORRFLUX cannot accept ',rec, ' as a parameter'
rec=-1
endif else rec=indgen(999) ;all records
endif
if n_elements(rec) eq 1 then rec=intarr(1)+rec ;convert to array
if rec(0) eq -1 then begin
helpme=1
goto,hlp
endif
;
setxy
if n_elements(dw) eq 0 then dw=100.
if w0def gt 0. then !x.range=[w0def-dw,w0def+dw]
pt='!6 CORRFLUX : '+file+':'
!x.title='!6Angstroms'
!y.title=ytit(0)
nrec=n_elements(rec)
for i=0,nrec-1 do begin
recno=rec(i)
gdat,file,h,w,f,e,recno
if w0def gt 0. then !x.range=[(w0def-dw)>w(0),(w0def+dw)<max(w)]
if n_elements(h) eq 1 then goto,done
!p.title=pt+strtrim(recno,2)+' '+strtrim(byte(h(100:139)>32b),2)
if ffl then begin
if n_elements(ddw) eq 0 then ddw=4
if ddw lt 0 then ddw=4
for ii=0,nl-1 do begin
k=fix(xindex(w,wl(ii))+0.5)
fl(ii)=total(f(k-ddw:k+ddw))/(w(k+ddw+1)-w(k-ddw))
endfor
cf=mean(fl)/fl ;correction factor
fact=interpol(cf,wl,w)
k1=fix(xindex(w,wl(0))+0.5)
k2=fix(xindex(w,wl(nl-1))+0.5)
fact(0:k1-ddw>0)=1.
fact(k2+ddw:*)=1.
if n_elements(nsm) eq 0 then nsm=5
if nsm le 0 then nsm=5
if nsm gt 2 then fact=smooth(fact,nsm)
endif else fact=interpol(fact0,w0,w)
fcor=f*fact
if not keyword_set(noplot) then begin
plot,w,f,psym=0
if (i eq 0) and (!d.name eq 'X') then wshow
oplot,w,fcor,psym=0,color=85
endif
if not keyword_set(NOSAVE) then kdat,file,h,w,fcor,e,rec(i) else print,' Data not saved'
endfor
done:
if keyword_set(stp) then stop,'CORRFLUX>>>'
return
end
corrflux.pro;2 000640 000451 000017 00000000000 06615503106 017056 1corrflux.pro ustar 00fwalter users 000002 260560 corrflux1.pro 000640 000451 000017 00000003256 06615503107 014422 0 ustar 00fwalter users 000002 260555 ;******************************************************************
function corrflux1,w,f
; called by CORRFLUX
if n_params(0) lt 2 then begin
print,' '
print,'* CORRFLUX1 is called by CORRFLUX - it is not a standalone procedure'
print,' '
return,0
endif
;
common comxy,xcur,ycur,zerr
plot,w,f
if !d.name eq 'X' then wshow
np=n_elements(f)-1
if (!x.range(0) eq 0) and (!x.range(1) eq 0) then ff=optfilt(f) else begin
k=where((w ge !x.crange(0)) and (w le !x.crange(1)))
k1=(k(0)-200)>0 & k2=(max(k)+200)<np
ff=f
z=optfilt(f(k1:k2))
ff(k1)=z
endelse
oplot,w,ff
xcur=mean(!x.crange)
ycur=mean(!y.crange)
if !d.name eq 'X' then wshow
print,string(7b)
print,' mark 2 extrema, Q to quit'
blowup,-1
x1=xcur
if (zerr eq 26) or (zerr eq 81) or (zerr eq 113) then return,-999
blowup,-1
x2=xcur
if (zerr eq 26) or (zerr eq 81) or (zerr eq 113) then return,-999
if x1 gt x2 then begin
t=x1
x1=x2
x2=t
endif
print,' mark region to exclude from interpolation, 0 if none'
blowup,-1
if zerr ne 48 then begin
x3=xcur
blowup,-1
x4=xcur
if x3 gt x4 then begin
t=x3
x3=x4
x4=t
endif
endif else begin
x3=-1 & x4=-1
endelse
k=where(w gt x1) & k1=k(0)>0
k=where(w gt x2) & k2=k(0)>0
k=where(w gt x3) & k3=k(0)>0
k=where(w gt x4) & k4=k(0)>0
length=k2-k1+1
val1=mean(ff((k1-10)>0:(k1+10)<np))
val2=mean(ff((k2-10)>0:(k2+10)<np))
v1=val1+(val2-val1)/length*findgen(length) ;linear interpolation
;
fact=ff
fact(k1)=v1
fact=fact/ff
if x3 ge 0 then begin ;avoidance region
val3=fact(k3)
val4=fact(k4)
length2=k4-k3+1
v2=val3+(val4-val3)/length2*findgen(length2) ;linear interpolation
fact(k3)=v2
endif
return,fact
end
corrflux1.pro;2 000640 000451 000017 00000000000 06615503107 017225 1corrflux1.pro ustar 00fwalter users 000002 260555 corrflux2.pro 000640 000451 000017 00000001226 06615503107 014420 0 ustar 00fwalter users 000002 260557 ;******************************************************************
function corrflux2,w,f
; called by FUDGEFLUX
if n_params(0) lt 2 then begin
print,' '
print,'* CORRFLUX2 is called by FUDGEFLUX - it is not a standalone procedure'
print,' '
return,0
endif
;
common comxy,xcur,ycur,zerr
plot,w,f
if !d.name eq 'X' then wshow
np=n_elements(f)-1
ff=optfilt(f,0,31,63,63,103)
oplot,w,ff
xcur=mean(!x.crange)
ycur=mean(!y.crange)
if !d.name eq 'X' then wshow
;
sm1=np/40 & if sm1 mod 2 eq 0 then sm1=sm1+1
sm2=np/20 & if sm2 mod 2 eq 0 then sm2=sm2+1
v1=smooth(maxfilt(f,-1,sm1),sm2)
oplot,w,v1,color=15
;
fact=v1/ff
print,string(7b)
return,fact
end
corrflux2.pro;2 000640 000451 000017 00000000000 06615503107 017231 1corrflux2.pro ustar 00fwalter users 000002 260557 createfile.pro 000640 000451 000017 00000000511 06615503107 014567 0 ustar 00fwalter users 000002 260564 ;*************************************************************************
pro createfile,name
if n_params(0) eq 0 then begin
print,' '
print,'* CREATEFILE - create an ascii file'
print,'* calling sequence: CREATEFILE,filename'
print,' '
return
endif
get_lun,lu
openw,lu,name
close,lu
free_lun,lu
return
end
createfile.pro;2 000640 000451 000017 00000000000 06615503107 017561 1createfile.pro ustar 00fwalter users 000002 260564 crosscor.pro 000640 000451 000017 00000005341 06615503107 014330 0 ustar 00fwalter users 000002 260565 ;**************************************************************************
pro crosscor,w1,f10,w2,f20,dw,xc,logl,cut,a,debug=stp,excl=badlam, $
NOSMOOTH=NOSMOOTH,FWID=FWID
common comxy,xcur,ycur,zerr,zzz,lu3
xcur=1. & ycur=0.
; cross correlate spectra f1(w1), f2(w2)
; return cross correlation product vector in XC
; if logl eq 0 (default), use log-lambda correlations
;
f1=f10 & f2=f20
if not keyword_set(stp) then stp=0
if not keyword_set(badlam) then badlam=-1
if n_elements(badlam) lt 2 then iexcl=0 else iexcl=1
if iexcl eq 1 then begin
nb=n_elements(badlam)
odd=nb mod 2
if odd eq 1 then badlam=badlam(0:nb-2)
endif
if n_params(0) lt 7 then logl=0 ;default to log-lambda scale
np=n_elements(f2)
if n_elements(cut) eq 0 then cut=(50<(np/3))
IF NOT KEYWORD_SET(NOSMOOTH) THEN BEGIN
fftsm,f1,1
fftsm,f2,1
ENDIF
if logl eq 0 then begin ;put on log-lambda scale
loglam,w1,f1,lw1,lf1
np=n_elements(lw1)
dw=(lw1(np-1)-lw1(0))/(np-1)
loglam,w2,f2,lw1,lf2,1
if iexcl eq 1 then bl=alog10(badlam)
endif else begin
lf1=f1 & lf2=f2
bl=badlam
endelse
;
if (n_elements(w1) le 1) or (n_elements(w2) le 1) then goto,nowave
if iexcl eq 1 then begin ;exclude data segments
nb=n_elements(bl)/2
if nb ge 1 then for i=0,nb-1 do begin
bw1=bl(i*2) & bw2=bl(i*2+1)
k=where((lw1 gt bw1) and (lw1 lt bw2))
if k(0) ne -1 then linterpl,lf1,(k(0)-1)>0,k(n_elements(k)-1)
if k(0) ne -1 then linterpl,lf2,(k(0)-1)>0,k(n_elements(k)-1)
endfor
endif
;
nowave:
IF N_ELEMENTS(FWID) EQ 0 THEN FWID=51
if fwid le 0 then fwid=51
filt1=smooth(maxfilt(lf1),FWID)
filt2=smooth(maxfilt(lf2),FWID)
ff1=lf1/filt1 ;divide out continuum
ff2=lf2/filt2
ff1=ff1-mean(ff1) ;subtract mean
ff2=ff2-mean(ff2)
if iexcl eq 1 then begin ;re-exclude data segments
if nb ge 1 then for i=0,nb-1 do begin
bw1=bl(i*2) & bw2=bl(i*2+1)
k=where((lw1 gt bw1) and (lw1 lt bw2))
if k(0) ne -1 then ff1(k)=0.
if k(0) ne -1 then ff2(k)=0.
endfor
endif
spxcor,ff1,ff2,xc,cut
;
if n_params(0) ge 9 then begin ;do not fit here
x=findgen(n_elements(xc))
xcc=xc
if stp eq 3 then begin
plot,x,xc
print,' put cursor on peak'
blowup,-1
c=xcur
xcc=xc
k=where(x lt c-15) & if k(0) ne -1 then xcc(k)=0
k=where(x gt c+15) & if k(0) ne -1 then xcc(k)=0
endif else begin
c=where(xc eq max(xc)) & c=c(0)
endelse
;
a=[max(xc((c-5)>0:(c+5)<(n_elements(x)-1)))-mean(xc),c,5.,mean(xc),0.,0.]
z=gaussfit(x,xcc,a)
a(1)=a(1)-max(x)/2. ;shift
endif
;
if stp gt 0 then begin
!p.position=[.2,.55,.9,.95]
;set_viewport,.2,.9,.55,.95
plot,lw1,ff1 & oplot,lw1,ff2,color=15
endif
return
end
crosscor.pro;2 000640 000451 000017 00000000000 06615503107 017046 1crosscor.pro ustar 00fwalter users 000002 260565 ctit.pro 000640 000451 000017 00000001300 06615503110 013413 0 ustar 00fwalter users 000002 260570 ;************************************************************************
pro ctit,h,nch,newt ;change title (words 100-159) in header
nh=n_elements(h)
nch=1 ; no change made
if nh lt 100 then return ;header too short
nh=(nh-100)<60
hmax=99+nh
if hmax lt 0 then return
title=string(byte(h(100:hmax)>32))
t=bytarr(nh)
if n_elements(newt) eq 0 then newt=''
if strlen(newt) eq 0 then begin
print,' old title:',title
newt=' '
read,' enter new title, CR if OK: ',newt
endif
if newt eq '' then return
if strlen(newt) gt nh then newt=strmid(newt,0,nh-1)
t(0)=byte(newt)
h(100)=t
nch=0
print,'revised title: ',string(byte(h(100:159)))
return
end
ctit.pro;2 000640 000451 000017 00000000000 06615503110 015250 1ctit.pro ustar 00fwalter users 000002 260570 cutgaps.pro 000640 000451 000017 00000000751 06615503110 014130 0 ustar 00fwalter users 000002 260571 ;*********************************************************************
pro cutgaps,h,w,f,e,ngap,wgap1,wgap2
h(900)=ngap
if ngap le 0 then return
help,w
for i=0,ngap-1 do begin
print,' cutting gap',i,wgap1(i),wgap2(i)
dw=wgap2(i)-wgap1(i)
k=where((w lt wgap1(i)) or (w gt wgap2(i)))
w=w(k)
f=f(k)
e=e(k)
h(901+i*4)=fix(wgap1(i))
h(902+i*4)=fix(30000.*(wgap1(i)-h(901+i*4)))
h(903+i*4)=fix(dw)
h(904+i*4)=fix(30000.*(dw-h(903+i*4)))
endfor
help,w
return
end
cutgaps.pro;2 000640 000451 000017 00000000000 06615503110 016457 1cutgaps.pro ustar 00fwalter users 000002 260571 cwhere.pro 000640 000451 000017 00000002443 06615503110 013740 0 ustar 00fwalter users 000002 260572 pro cwhere,w,accurate=accurate,pixels=pixels
common comxy,xcur,ycur,zerr,RESETSCALE,lu3
common icurunits,xunits,yunits,title,c1,c2,c3,ch,c4,c5,c6,c7,c8,c9
if strupcase(!d.name) ne 'X' then return
wshow
ylog=!y.type
xlog=!x.type
yl=!y.crange(0) & yu=!y.crange(1)
if ylog eq 1 then begin
yl=10^yl & yu=10^yu
endif
xl=!x.crange(0) & xu=!x.crange(1)
if xlog eq 1 then begin
xl=10^xl & xu=10^xu
endif
if n_elements(w) lt 2 then pixels=0
;
print,' Use left mouse button to stop, center to mark, right to mark and print'
cr=string("15b)
if n_elements(xunits) eq 1 then sx=xunits+': ' else sx='X: '
if n_elements(yunits) eq 1 then sy=yunits+': ' else sy='Y: '
form="($,A11,F9.3,', ',A9,G11.3,a)"
if keyword_set(accurate) then form="($,A11,F12.6,', ',A9,G14.6,a)"
zerr=0
print,format="($,A)",string("12b) ;lf
while zerr lt 1 do begin
cursor,x,y,2
zerr=!err
if keyword_set(pixels) then v=xindex(w,x) else v=x
if ((x ge xl) and (x le xu) and $
(y ge yl) and (y le yu)) then print,form=form,sx,v,sy,y,cr
if zerr eq 4 then begin
print,format="($,A)",string("12b)
zerr=0
tkp,1,x,y
wait,0.5
endif
if zerr eq 2 then begin
zerr=0
tkp,1,x,y
wait,0.5
endif
endwhile
print,format="($,A)",string("12b)
xcur=x
ycur=y
zerr=119
return
end
cwhere.pro;2 000640 000451 000017 00000000000 06615503110 016076 1cwhere.pro ustar 00fwalter users 000002 260572 dat_to_icd.pro 000640 000451 000017 00000003263 06615503110 014557 0 ustar 00fwalter users 000002 260574 ;******************************************************************************
pro dat_to_icd,in,recs,out=out,maxrec=maxrec,helpme=helpme,notitle=notitle
if n_params(0) le 0 then helpme=1
if keyword_set(helpme) then begin
print,' '
print,'* DAT_TO_ICD'
print,'* convert old format ICUR data files to new format files'
print,'* '
print,'* calling sequence:
print,'* DAT_TO_ICD,infile,rec1,rec2'
print,'* INFILE: name of input .DAT data file'
print,'* RECS: optional list of records to be copied, default=all'
print,'*'
print,'* KEYWORDS:'
print,'* OUT: name of output file, default=INFILE.ICD'
print,'* MAXREC: set to force vector length to max in input file'
print,'* NOTITLE: set to override query for title'
print,' '
return
endif
;
infile=in
outfile=in
if get_ext(in) eq '' then infile=infile+'.dat' else begin
k=strpos(in,'.')
outfile=strmid(in,0,k)
endelse
if not ffile(infile) then begin
bell,3
print,' file ', infile,' not found'
return
endif
;
if keyword_set(out) then outfile=out
if keyword_set(notitle) then qt=1 else qt=0
;
if n_params(0) lt 2 then recs=indgen(999)
nr=n_elements(recs)
if keyword_set(maxrec) then begin
reclen=0
for i=0,nr-1 do begin
getdat,in,h,w,f,e,recs(i)
if n_elements(h) lt 2 then goto,done1
reclen=h(7)>reclen
endfor
done1:
print,'reclen=',reclen
endif
;
for i=0,nr-1 do begin
getdat,in,h,w,f,e,recs(i)
if n_elements(h) lt 2 then goto,done
if keyword_set(maxrec) then kdat,outfile,h,w,f,e,vlen=reclen,notitle=qt $
else kdat,outfile,h,w,f,e,notitle=qt
endfor
done:
ldat,outfile
print,' New format data file is ',outfile
return
end
dat_to_icd.pro;2 000640 000451 000017 00000000000 06615503110 017530 1dat_to_icd.pro ustar 00fwalter users 000002 260574 deconv_ghrs.pro 000640 000451 000017 00000004077 06615503111 014762 0 ustar 00fwalter users 000002 260600 ;**********************************************************************
function deconv_ghrs,w,f,grat,psf_file=psf_file
common com1,h
COMMON COMXY,XCUR,YXUR,ZERR
common icurunits,xunits,yunits,title,c1,c2,c3,ch
; psfs for 120,140,170,210,270
c=[120,140,170,210] ;,270] 270 not available
;
if not keyword_set(psf_file) then begin
wcen=mean(w)
d=abs(wcen-c*10.)
m=where(d eq min(d))
cent=c(m) & CENT=CENT(0)
np=n_elements(f)
case 1 of
np le 600: ext='_1'
np ge 1500: ext='_4'
else: ext='_2'
endcase
file='psf'+strtrim(cent,2)+ext
endif else file=psf_file
;
icurdata=getenv('icurdata')
file=icurdata+file
if get_ext(file) eq '' then file=file+'.dat'
print,' PSF FILE: ',file
s=0. & psf=0.
openr,lu,file,/get_lun
genrd,lu,s,psf
;psf0=psf ;units are arcsec
;
; disp is dispersion in A/diode
dgrat=[.054,.069,.078,.092,.57] ;crude dispersions
if (grat ge 1) and (grat le 5) then begin
disp=dgrat(grat-1)
dpa=disp*4. ;angstroms per arcsec ; diode=0.25 arcsec
s=s*dpa ;PSF scale per angstrom
ds=max(s)-min(s) ;length of psf
dl=w(1)-w(0) ;bin size
npsfp=fix(ds/dl) ;number of points
if npsfp mod 2 eq 0 then npsfp=npsfp+1
hp=(npsfp-1)/2
ss=(findgen(npsfp)-hp)*dl
psf=interpol(psf,s,ss)
endif
; disp=0.
; print,' WARNING: Dispersion undefined. please enter dispersion in A/diode'
; read,disp
; endelse
;
psf=psf/total(psf)
print,' type number of iterations (1-9, 0=10), Q to abort'
blowup,-1
if (zerr lt 48) or (zerr gt 57) then return,f
niter=zerr-48
if niter eq 0 then niter=10
lucy_guess,niter,f,psf,newf
h(2)=1000+niter
;
while (niter ge 0) and (niter le 10) do begin
oplot,w,newf,color=c2
nd=strtrim((h(2)-1000),2)
print,nd,' iterations done. Enter number to do (1-9, 0=10), Q to abort'
blowup,-1
if (zerr lt 26) or (zerr gt 81) or (zerr eq 113) then return,f
if (zerr lt 48) or (zerr gt 57) then return,newf
niter=zerr-48
if niter eq 0 then niter=10
lucy_guess,niter,f,psf,newf,guess=newf
h(2)=h(2)+niter
endwhile
;
return,newf
end
deconv_ghrs.pro;2 000640 000451 000017 00000000000 06615503111 020137 1deconv_ghrs.pro ustar 00fwalter users 000002 260600 degrade.pro 000640 000451 000017 00000002243 06615503111 014061 0 ustar 00fwalter users 000002 260576 ;************************************************************************
PRO DEGRADE,W,F,E,BW,BF ; DEGRADE IUE HI-RES DATA INTO LO-RES DATA
; W,F,E ARE OVERWRITTEN, AND IDAT IS SET TO 0 (LO-RES)
COMMON COM1,H,IK,IFT,NSM,C
common comxy,xcur,ycur,zerr
csave=c & nsave=nsm
IF h(3) gt 4 THEN BEGIN
PRINT,' You cannot DEGRADE this type of data.'
zerr=68 ;do plot
RETURN
ENDIF
S=N_ELEMENTS(W)
T1=F
T2=W
DISP=(W(50)-W(0))/50.
if disp gt 1. then begin
print,' This is not high resolution data.'
return
endif
case 1 of
h(3) ge 3: lodisp=1.2 ;SW
(h(3) eq 1) or (h(3) eq 2): lodisp=1.88 ;LW
else: begin
print,'DEGRADE: invalid camera - returning'
return
end
endcase
IDIM=FIX(lodisp/DISP)
if idim le 1 then begin
print,'DEGRADE: dispersion ratios =',idim,' Returning'
return
endif
;
V=8.0 ;LWP
IF h(3) ge 3 THEN V=6.0 ;SWP
ROTVEL,-2,W,V/2.
S=N_ELEMENTS(C)
IF S gt 1 THEN F=CONVOL(F,C)
c=csave & nsm=nsave
coadd,w,f,e,idim
;
H(51)=FIX((V-10000)/10.)
H(52)=100*FIX(V-FIX(V))
H(1)=2
H(2)=999
BDATA,H,-1,W,F,E,BW,BF
ZERR=68
RETURN
END
degrade.pro;2 000640 000451 000017 00000000000 06615503111 016337 1degrade.pro ustar 00fwalter users 000002 260576 degtodms.pro 000640 000451 000017 00000001715 06615503111 014266 0 ustar 00fwalter users 000002 260603 ;*************************************************************************
pro DEGtoDMS,d,dd,dm,ds,prt=prt,nsigfig=nsigfig
s=size(d)
ndim=s(0)
if ndim eq 0 then d=intarr(1)+d
np=n_elements(d)
dd=fix(d)
dm=fix((d-dd)*60)
ds=(d-dd-dm/60.)*3600.
k=where(d lt 0)
if k(0) ne -1 then begin ;negative declinations - leave only most sign. neg
k1=where(dd lt 0)
if k1(0) ne -1 then begin
dm(k1)=abs(dm(k1)) & ds(k1)=abs(ds(k1))
endif
k1=where(dm lt 0)
if k1(0) ne -1 then ds(k1)=abs(ds(k1))
endif
if ndim eq 0 then begin
d=d(0)
dd=dd(0)
dm=dm(0)
ds=ds(0)
endif
if n_params(0) eq 1 then prt=1
hfmt='(I4)'
mfmt='(I3)'
if n_elements(nsigfig) eq 0 then nsigfig=2
nf=nsigfig+4
fmt='(F'+strtrim(nf,2)+'.'+strtrim(nsigfig,2)+')'
if keyword_set(prt) then begin
if ndim eq 0 then print,string(dd,hfmt),string(dm,mfmt),string(ds,fmt) $
else for i=0,np-1 do print,string(dd(i),hfmt),string(dm(i),mfmt),string(ds(i),fmt)
endif
return
end
degtodms.pro;2 000640 000451 000017 00000000000 06615503111 016754 1degtodms.pro ustar 00fwalter users 000002 260603 degtohms.pro 000640 000451 000017 00000001161 06615503112 014267 0 ustar 00fwalter users 000002 260604 ;*************************************************************************
pro DEGtoHMS,a,h,m,s,prt=prt,nsigfig=nsigfig
s=size(a)
ndim=s(0)
np=n_elements(a)
k=where(a lt 0.,nk) & if nk gt 0 then a(k)=360.+a(k)
h=fix(a/15.)
m=fix((a-h*15)*4.)
s=(a-(h*15.)-m/4.)*240.
if n_params(0) eq 1 then prt=1
if n_elements(nsigfig) eq 0 then nsigfig=3
nf=4+nsigfig
fmt='(F'+strtrim(nf,2)+'.'+strtrim(nsigfig,2)+')'
ifmt='(I3)'
if keyword_set(prt) then begin
if ndim eq 0 then print,string(h,ifmt),string(m,ifmt),string(s,fmt) else $
for i=0,np-1 do print,string(h(i),ifmt),string(m(i),ifmt),string(s(i),fmt)
endif
return
end
degtohms.pro;2 000640 000451 000017 00000000000 06615503112 016766 1degtohms.pro ustar 00fwalter users 000002 260604 deltamag.pro 000640 000451 000017 00000010262 06615503112 014232 0 ustar 00fwalter users 000002 260610 ;*******************************************************************
PRO deltamag,file,recs,w0,dw1,dw2,prt=prt
if n_params(0) eq 0 then file=-1
if ifstring(s) ne 1 then begin ;integer file passed
if file eq -1 then begin ;help
print,' '
print,'* DELTAMAG - narrow band magnitude index
print,'* calling sequence: DELTAMAG,file,recs,w0,dw1,dw2,'
print,'* file: name of file, no default'
print,'* recs: list or records to measure, default(-1)=all in file'
print,'* w0: central wavelength, defaults: (0) = 6563.'
print,'* w0=-1 for Stocke''s 6495A correlation'
print,'* dw1: narrow band width, default (0) = 5A'
print,'* dw2: broad band width, default (0) = 50A'
print,'*'
print,'* KEYWORDS:
print,'* IPRT: 1 to print listing; 0 (default) for noprint'
print,'* 2 generates deltamag.dta file (input to deltamagcor)'
print,'* output is in file DELTAMAG.LST on disk'
print,'*'
print,'* magnitude index = f(narrow)/f(broad)
print,' '
return
endif
endif
k=strpos(string(file),'jacobyatlas')
if k ne -1 then ibmv=1 else ibmv=0 ;1 if B-V stored in header
if n_params(0) lt 2 then recs=indgen(999)
if n_elements(recs) eq 1 then recs=intarr(1)+recs
if recs(0) eq -1 then recs=indgen(999)
;
if n_params(0) lt 3 then w0=0
case 1 of
w0 eq 0: begin
w0=6563. & dw1=5. & dw2=50.
idef=0
end
w0 eq -1: begin
w0=6495. & dw1=5. & dw2=50.
idef=1
end
else: begin
idef=-999
if n_params(0) lt 4 then dw1=5.
if dw1 lt 0. then dw1=5.
if n_params(0) lt 5 then dw2=50.
if dw2 lt 0. then dw2=50.
end
endcase
;
if dw1 gt dw2 then begin ;flip band widths
t=dw1
dw1=dw2
dw2=t
endif
if dw1 eq dw2 then begin
print,' band widths equal - null measurement'
return
endif
;
recs=fix(recs)
n=n_elements(recs)
if n eq 999 then print,' Measuring all records' else PRINT,N,' records to measure'
L=[w0-dw1,w0+dw1,w0-dw2,w0+dw2]
get_lun,lu
k=strpos(file,']')
if k eq -1 then ff=file else ff=strmid(file,k+1,32)
k=strpos(ff,'.')
if k ne -1 then ff=strmid(ff,0,k)
outfile='dm'+ff+strtrim(string(fix(w0)),2)+'.lst'
OPENW,lu,outfile
printf,lu,' Central wavelength:',w0,'A'
if idef eq 0 then printf,lu,' H-alpha dM/log B-V calibration'
if idef eq 1 then printf,lu,' John Stocke''s lambda 6495/B-V calibration'
PRINTF,lu,' Narrow band: ',L(0),L(1),' width=',dw1,'A'
PRINTF,lu,' Broad band: ',L(2),L(3),' width=',dw2,'A'
if n_params(0) gt 1 then printf,lu,' Data from file ',file
if iprt eq 2 then begin
OPENW,lu2,'DELTAMAG.dta',get_lun
printf,lu2,w0,dw1,dw2
endif
printf,lu,' '
FOR I=0,N-1 DO BEGIN
GDAT,file,H,W,F,E,recs(i)
if n_elements(h) eq 1 then goto,done ;no more records
LAB=STRTRIM(BYTE(H(100:139)>32),2)
NW=N_ELEMENTS(W)
IF (L(2) LT W(0)) OR (L(3) GT W(NW-1)) THEN goto,skip
indx=fix(xindex(w,l)+0.5) ;TABINV,W,L,INDX
IN=INDX(1)-INDX(0)
DN=W(INDX(1))-W(INDX(0))
FN=TOTAL(F,INDX(0),IN)/DN
IB=INDX(3)-INDX(2)
DB=W(INDX(3))-W(INDX(2))
FB=TOTAL(F,INDX(2),IB)/DB
DM=-2.5*ALOG10(FN/FB)
bmv0=float(h(63))/1000.
r=string(recs(i),'(I3)')
case 1 of
idef eq 0: begin ;Halpha
bmv=10^(0.18-2.17*dm)
dbmv=10^(0.006+0.06*dm)
dm=string(dm,'(F4.2)')
bmv=string(bmv,'(F5.2)')
printf,lu,' Rec:',r,' dM = ',dm,' B-V=',bmv,' +/- ',dbmv,' ',lab
end
idef eq 1: begin ; 6495
bmv=0.358+6.415*dm
dbmv=0.041+0.47*dm
dm=string(dm,'(F4.2)')
bmv=string(bmv,'(F5.2)')
dbmv='0.15'
printf,lu,' Rec:',r,' dM = ',dm,' B-V=',bmv,' +/- ',dbmv,' ',lab
end
else: begin
PRINTF,lu,' Rec:',R,' delta magnitude (narrow-broad) =',DM,' ',LAB
end
endcase
if iprt eq 2 then begin
if ibmv eq 1 then printf,lu2,recs(i),' ',dm,bmv0,' ',lab $
else printf,lu2,recs(i),' ',dm,' ',lab
endif
SKIP:
ENDFOR
done: CLOSE,lu
free_lun,lu
if iprt eq 2 then begin
close,lu2
free_lun,lu2
endif
print,' Output file is named ',outfile
if abs(iprt) eq 1 then spawn_print,outfile
RETURN
END
deltamag.pro;2 000640 000451 000017 00000000000 06615503112 016673 1deltamag.pro ustar 00fwalter users 000002 260610 dmstodeg.pro 000640 000451 000017 00000001327 06615503112 014270 0 ustar 00fwalter users 000002 260605 ;**********************************************************************
function dmstodeg,d,dm,ds
if n_params(0) eq 2 then ds=0.
if n_params(0) eq 1 then begin
if n_elements(d) lt 3 then return,-999
dm=abs(d(1)) & ds=abs(d(2)) & d=d(0)
endif
s=size(d)
if s(0) eq 0 then begin ;scalar
IF D NE 0 THEN ISIGN=D/ABS(D) ELSE BEGIN ;DEC=+/-0
IF DM NE 0 THEN ISIGN=DM/ABS(DM) ELSE BEGIN ;MIN=0
IF DS NE 0 THEN ISIGN=DS/ABS(DS) ELSE ISIGN=1.
ENDELSE
ENDELSE
endif else begin ;array
isign=fltarr(n_elements(d))+1.
k=where((d lt 0) or (dm lt 0) or (ds lt 0.))
if k(0) ne -1 then isign(k)=-1.
endelse
DELTA=ISIGN*(ABS(D)+ABS(DM)/60.+ABS(DS)/3600.)
return,delta
end
dmstodeg.pro;2 000640 000451 000017 00000000000 06615503112 016757 1dmstodeg.pro ustar 00fwalter users 000002 260605 dnd.pro 000640 000451 000017 00000001032 06615503112 013221 0 ustar 00fwalter users 000002 260606 ;***************************************************************
PRO DND,H,WAVE,FLUX ; DIVIDE BY ND RESPONSE
common comxy,xcur,ycur,zerr
common icdisk,icurdisk,icurdata
IF (h(3) lt 4) or (h(3) ge 100) THEN begin ;only valid for KPNO blue data
zerr=68
return ;wrong data type
endif
gdat,icurdata+'optstd',hs,w,f,e,0
IF ABS(W(0)-WAVE(0)) GT 500. THEN BEGIN
PRINT,' *** WAVELENGTH RANGE WRONG FOR QUARTZ REMOVAL'
zerr=68
RETURN
ENDIF
FF=INTERPOL(F,W,WAVE)
FLUX=FLUX*FF
RETURN
END
dnd.pro;2 000640 000451 000017 00000000000 06615503112 014656 1dnd.pro ustar 00fwalter users 000002 260606 dqcheck.pro 000640 000451 000017 00000001750 06615503113 014067 0 ustar 00fwalter users 000002 260607 ;************************************************************************
function dqcheck,fin,eps,q ;check data quality and censor at level q
if n_params(0) lt 2 then fin=-1
if n_elements(fin) le 1 then begin
print,' '
print,'DQcheck - check data quality vector and interpolate over bad points'
print,' calling sequence: FOUT=dqcheck(fin,eps,q)'
print,' FOUT: Output vector, interpolated over flagged points'
print,' FIN: Input vector'
print,' EPS: Data quality vector'
print,' Q: threshold. Points with EPS>Q are interpolated, def=0'
print,' '
return,fin
endif
if n_params(0) lt 3 then q=0
k=where(eps ge q,nk)
if nk eq 0 then begin
print,'DQCHECK: warning - no points with acceptable data quality'
print,' returning input vector'
return,fin ;no bad data
endif
np=n_elements(fin)
k=where(eps le q,nk)
if nk eq -1 then return,fin
x=findgen(np)
xx=x(k)
ff=fin(k)
fout=interpol(ff,xx,x)
return,fout
end
dqcheck.pro;2 000640 000451 000017 00000000000 06615503113 016352 1dqcheck.pro ustar 00fwalter users 000002 260607 drlin.pro 000640 000451 000017 00000002604 06615503113 013573 0 ustar 00fwalter users 000002 260615 ;*****************************************************
pro drlin,w,ls=ls,COLOR=COLOR,vertical=vertical,helpme=helpme
; PROCEDURE DRLIN TO DRAW LINE AT LEVEL W
if keyword_set(helpme) then begin
print,' '
print,'* DRLIN - draw horizontal line on plot '
print,'* calling sequence: DRLIN,Y '
print,'* Y: height of line, default=0'
print,'* '
print,'* KEYWORD:'
print,'* LS: linestyle value, default=0'
print,'* COLOR: color of line, default=!p.color'
print,'* VERTICAL: draw vertical line'
print,' '
return
endif
if n_elements(w) eq 0 then w=0. ;default
w=w(0)
;
if keyword_set(vertical) then begin
yt=!x.type
yrange=!x.crange
xt=!y.type
xrange=!y.crange
endif else begin
yt=!y.type
xt=!x.type
yrange=!y.crange
xrange=!x.crange
endelse
;
if (yt eq 1) and (w le 0.) then return ; out of bounds
y1=yrange(0) & y2=yrange(1)
if y2 lt y1 then begin
y2=yrange(0) & y1=yrange(1)
endif
case 1 of ;return if out of bounds
yt eq 1: IF (alog10(W) GT Y2) OR (alog10(W) LT Y1) THEN RETURN
else: IF (W GT Y2) OR (W LT Y1) THEN RETURN
endcase
if xt eq 1 then x=10^(xrange) else X=xrange
Y=X*0.+W
if not keyword_set(ls) then ls=0
if not keyword_set(color) then color=!p.color
if keyword_set(vertical) then begin
t=x & x=y & y=t
endif
OPLOT,X,Y,psym=0,linestyle=ls,color=color
return
END
drlin.pro;2 000640 000451 000017 00000000000 06615503113 015565 1drlin.pro ustar 00fwalter users 000002 260615 ebyte.pro 000640 000451 000017 00000001255 06615503113 013576 0 ustar 00fwalter users 000002 260617 ;****************************************************************
function ebyte,e
n=n_elements(e)
if n gt 32767 then e1=bytarr(32768L) else e1=bytarr(n)
k=where(e gt 0) & if k(0) ne -1 then e1(k)=0b ;good data
k=where((e gt -201) and (e le 0)) & if k(0) ne -1 then e1(k)=1b ;extrapolated ITF
k=where((e gt -221) and (e le -201)) & if k(0) ne -1 then e1(k)=2b ;microphonics
k=where((e gt -301) and (e le -221)) & if k(0) ne -1 then e1(k)=3b ;hot spot
k=where((e gt -801) and (e le -301)) & if k(0) ne -1 then e1(k)=4b ;reseau
k=where((e gt -1601) and (e le -801)) & if k(0) ne -1 then e1(k)=5b ;saturated
k=where(e le -1601) & if k(0) ne -1 then e1(k)=6b ;uncorrected
return,e1
end
ebyte.pro;2 000640 000451 000017 00000000000 06615503113 015567 1ebyte.pro ustar 00fwalter users 000002 260617 eqwid.pro 000640 000451 000017 00000004566 06615503113 013601 0 ustar 00fwalter users 000002 260620 ;**************************************************************************
function eqwid,w,f,lamcen,dl,b1,b2,db1,db2,stp=stp ;estimate equivalent widths
if n_params(0) lt 2 then begin
print,' '
print,'* function EQWID - estimate equivalent widths'
print,'* calling sequence X=EQWID(W,F,LAM,DL,B1,B2,DB1,DB2)'
print,'* W,F: input wavelength and flux vectors'
print,'* LAM: wavelength of feature to be measured, def = H-alpha'
print,'* DL: width to integrate over, default=5 A'
print,'* B1,B2: beginning wavelengths of background region,'
print,'* defaults=LAM-20 A,LAM+10 A'
print,'* DB1,BD2: widths of background region, default=10 Angstroms'
print,'* X: output equivalent width in units of W'
print,' '
return,0.
endif
dw=mean(w(1:*)-w) ;Angstroms per bin
;
if n_params(0) lt 3 then lamcen=-1
if n_params(0) lt 4 then dl=-1
if n_params(0) lt 8 then db2=-1
if n_params(0) lt 7 then db1=-1
if n_params(0) lt 6 then b2=-1
if n_params(0) lt 5 then b1=-1
;
if lamcen lt 0 then lamcen=6563. ;default to H-alpha
if dl lt 0 then dl=5. ;default to 5 Angstroms width
if db1 lt 0 then db1=10. ;default to 20 bins
if db2 lt 0 then db2=db1 ;default to 20 bins
if b2 lt 0 then b2=lamcen+10. ;default to Lamcen+20 bins
if b1 lt 0 then b1=lamcen-20. ;default to Lamcen-10 bins
if (b1 lt w(0)) or (b2+db2 gt max(w)) then begin
print,' region not completely contained within spectrum - returning'
return,-9999.
endif
;
ff=optfilt(f)
dl2=dl/2.
WV=[b1,b1+db1,b2,b2+db2,lamcen-dl2,lamcen+dl2]
ii=xindex(w,wv) ;tabinv,w,wv,ii
; do first background region
mw1=b1+db1/2. ;mean wavelength
fi1=fix(ii(0)) & fi2=fix(ii(1))
fb1=total(ff(fi1+1:fi2))+ff(fi1)*(fi1+1.-ii(0))+ff(fi2+1)*(ii(1)-fi2)
fb1=fb1*dw/(db1+dw) ;flux per A
; do second background region
mw2=b2+db2/2. ;mean wavelength
fi1=fix(ii(2)) & fi2=fix(ii(3))
fb2=total(ff(fi1+1:fi2))+ff(fi1)*(fi1+1.-ii(2))+ff(fi2+1)*(ii(3)-fi2)
fb2=fb2*dw/(db2+dw) ;flux per A
;
; do line
fi1=fix(ii(4)) & fi2=fix(ii(5))
fl=total(f(fi1+1:fi2))+f(fi1)*(fi1+1.-ii(4))+f(fi2+1)*(ii(5)-fi2)
fl=fl*dw
;
fb=fb1+(fb2-fb1)*(lamcen-mw1)/(mw2-mw1) ;interpolated background
ew=(fb*dl-fl)/fb
;
if keyword_set(stp) then stop,'EQWID>>>'
return,ew
end
eqwid.pro;2 000640 000451 000017 00000000000 06615503113 015563 1eqwid.pro ustar 00fwalter users 000002 260620 erbar.pro 000640 000451 000017 00000005173 06615503114 013562 0 ustar 00fwalter users 000002 260623 ;*************************************************************************
PRO ERBAR,IXY,a1,a2,a3,a4,a5,a6,color=color,helpme=helpme
; ixy=1 (X), 2 (Y), 3 (X+Y)
;
if n_params(0) lt 3 then helpme=1
if keyword_set(helpme) then begin
print,' '
print,' * ERBAR - overplot error bars'
print,' * calling sequences: ERBAR,i,y,y1,y2,x,x1,x2'
print,' * ERBAR,i,x,dx,y,dy'
print,' * ERBAR,i,z,z1,z2,w'
print,' * ERBAR,i,z,dz,w'
print,' * I: 1: error bar in X'
print,' * I: 2: error bar in Y'
print,' * I: 3: error bars in both coordinates'
print,' * X,Y: X,Y points'
print,' * dX,dY : 1 sigma error'
print,' * y1,y2,x1,x2: actual values of ends of error bars'
print,' * Z,W: X or Y, depending on I
print,' *'
print,' * KEYWORDS:'
print,' * COLOR: plot color'
print,' '
return
endif
if (ixy le 0) or (ixy gt 3) then return
case 1 of
(n_params(0) eq 4) or (n_params(0) eq 3): begin
if ixy eq 3 then return ;invalid
if ixy eq 1 then begin
x0=a1
x1=a1-a2 & x2=a1+a2
if n_params(0) eq 4 then y0=a3 else y0=indgen(n_elements(x0))
y1=y0 & y2=y0
endif else begin
y0=a1 & y1=a1-a2 & y2=a1+a2
if n_params(0) eq 4 then x0=a3 else x0=indgen(n_elements(y0))
x1=x0 & x2=x0
endelse
end
(n_params(0) eq 5) and (ixy lt 3): begin
if ixy eq 1 then begin
x0=a1 & x1=a2 & x2=a3
y0=a4 & y1=a4 & y2=a4
endif else begin
y0=a1 & y1=a2 & y2=a3
x0=a4 & x1=a4 & x2=a4
endelse
end
(n_params(0) eq 5) and (ixy eq 3): begin
x0=a1 & x1=a1-a2 & x2=a1+a2
y0=a3 & y1=a3-a4 & y2=a3+a4
end
n_params(0) eq 7 : begin
ixy=3
x0=a1 & x1=a2 & x2=a3
y0=a4 & y1=a5 & y2=a6
end
else: begin
print,' Invalid input'
return
end
endcase
;
ix=ixy and 1
iy=ixy and 2
if !y.type eq 1 then y1=y1>10^(!y.crange(0))
;
NST=N_ELEMENTS(X0)
ny=n_elements(y)
if not keyword_set(color) then color=!p.color
if nst gt ny then begin
if n_elements(y0) gt 0 then y0=[y0,replicate(0,nst-ny)]
if n_elements(y1) gt 0 then y1=[y1,replicate(0,nst-ny)]
if n_elements(y2) gt 0 then y2=[y2,replicate(0,nst-ny)]
endif
FOR I=0L,NST-1L DO BEGIN
if iy gt 0 then begin
x=[x0(i),x0(i)]
y=[y1(i),y2(i)]
oplot,x,y,psym=0,color=color
endif
if ix gt 0 then begin
x=[x1(i),x2(i)]
; x=[x1(i)>!x.crange(0),x2(i)<!x.crange(1)]
y=[y0(i),y0(i)]
oplot,x,y,psym=0,color=color
endif
ENDFOR
RETURN
END
erbar.pro;2 000640 000451 000017 00000000000 06615503114 015533 1erbar.pro ustar 00fwalter users 000002 260623 errdiv.pro 000640 000451 000017 00000002653 06615503114 013765 0 ustar 00fwalter users 000002 260626 ;**************************************************************************
function ERRDIV,X0,DX0,Y0,DY0,print=print,unc=unc
;
if n_params(0) lt 2 then begin
print,' '
print,' * ERRDIV -- return statistical error on X/Y'
print,' * calling sequence: UNC=ERRDIV(X,dX,Y,dY)'
print,' * X: numerator'
print,' * dX: uncertainty on numerator'
print,' * Y: denominator'
print,' * dY: uncertainty on denominator'
print,' * KEYWORDS:'
print,' * PRINT: set to print result to screen'
print,' * UNC: set for relkative uncertainty. E.g., UNC=0.1 for 10% error'
print,' *'
print,' * call as UNC=ERRDIV(X,Y) to use sqrt(X),sqrt(Y) as errors'
print,' '
return,-999
endif
;
x=x0 & dx=dx0
fx=float(x)
if n_elements(unc) eq 0 then unc=0
if n_elements(unc) gt 0 then unc=unc(0)
unc=abs(unc)
case 1 of
unc gt 1.: begin
if unc lt 100. then unc=unc/100. else unc=0 ;change percent to fraction
end
else:
endcase
case 1 of
unc gt 0: begin
fy=float(dx)
dx=fx*unc
dy=fy*unc
end
n_params(0) eq 2: begin
fy=float(dx)
dx=sqrt(fx)
dy=sqrt(fy)
end
n_params(0) eq 3: begin
y=y0
fy=float(y)
dy=sqrt(fy)
end
else: begin
y=y0
dy=dy0
fy=float(y)
endelse
endcase
;
err=SQRT(DX*DX+(fX*DY/fY)*(fX*DY/fY))/fY
if keyword_set(print) then print,fx/fy,'+\-',err
return,err
end
errdiv.pro;2 000640 000451 000017 00000000000 06615503114 016136 1errdiv.pro ustar 00fwalter users 000002 260626 euler.pro 000640 000451 000017 00000005044 06615503114 013576 0 ustar 00fwalter users 000002 260630 ;************************************************************************
PRO EULER,AI,BI,AO,BO,SELECT ; from IDL USERs S/W
;+
; NAME:
; EULER
; PURPOSE:
; Transform between galactic, celestial, and ecliptic coordinates.
; Use the procedure ASTRO to use this routine interactively
; CALLING SEQUENCE:
; EULER,AI,BI,AO,BO,[SELECT]
; INPUTS:
; AI - Input Longitude in DEGREES, scalar or vector. If only two parameters
; are supplied, then AI and BI will be modified to contain the output
; longitude and latitude.
; BI - Input Latitude in DEGREES
; OPTIONAL INPUT:
; SELECT - Integer (1-6) specifying type of coordinate transformation.
;
; SELECT From To | SELECT From To
; 1 RA-DEC(1950) GAL.(ii) | 4 ECLIPTIC RA-DEC
; 2 GAL.(ii) RA-DEC | 5 ECLIPTIC GAL.(ii)
; 3 RA-DEC ECLIPTIC | 6 GAL.(ii) ECLIPTIC
;
; If omitted, program will prompt for the value of SELECT
; OUTPUTS:
; AO - Output Longitude in DEGREES
; BO - Output Latitude in DEGREES
; REVISION HISTORY:
; Written W. Landsman, February 1987
; Adapted from Fortran by Daryl Yentis NRL
;-
twopi = 6.2831853072d
fourpi = 12.566370614d
cdr = 360./twopi
;
psi = [ 0.57595865315D, 4.9261918136D, $
0.00000000000D, 0.0000000000D, $
0.11129056012D, 4.7005372834D]
stheta =[ 0.88781538514D,-0.88781538514D, $
0.39788119938D,-0.39788119938D, $
0.86766174755D,-0.86766174755D]
ctheta =[ 0.46019978478D, 0.46019978478D, $
0.91743694670D, 0.91743694670D, $
0.49715499774D, 0.49715499774D]
phi = [ 4.9261918136D, 0.57595865315D, $
0.0000000000D, 0.00000000000D, $
4.7005372834d, 0.11129056012d]
;
npar = n_params(0)
if npar lt 5 then begin
print,' '
print,' 1 RA-DEC(1950) TO GAL.(ii)
print,' 2 GAL.(ii) TO RA-DEC
print,' 3 RA-DEC TO ECLIPTIC
print,' 4 ECLIPTIC TO RA-DEC
print,' 5 ECLIPTIC TO GAL.(ii)
print,' 6 GAL.(ii) TO ECLIPTIC
;
read,'Enter selection',select
endif
I = select - 1 ; IDL offset
a = ai/cdr - phi(i)
b = bi/cdr
sb = sin(b) & cb = cos(b)
cbsa = cb * sin(a)
b = -stheta(i) * cbsa + ctheta(i) * sb
bo = asin(b<1.0d)*cdr
;
a = atan(ctheta(i) * cbsa + stheta(i) * sb, cb * cos(a))
;
; factor of 1./(cos(bo)) removed from both sin(a) and cos(a)
;
ao = ((a+psi(i)+fourpi) mod twopi) * cdr
if npar eq 2 then begin
ai = ao & bi=bo
endif
return
end
euler.pro;2 000640 000451 000017 00000000000 06615503114 015573 1euler.pro ustar 00fwalter users 000002 260630 euvespec.pro 000640 000451 000017 00000006440 06615503115 014306 0 ustar 00fwalter users 000002 260633 ;****************************************************************************
pro euvespec,lam,d,icd=icd,helpme=helpme,stp=stp,sw=sw,mw=mw,lw=lw, $
allspec=allspec,title=title,file=file
if keyword_set(helpme) then begin
print,' '
print,'* EUVESPEC - transform EUVE *_SPEC.FITS files into .ICD format'
print,'* calling sequence: EUVESPEC,W,F'
print,'* W,F: output wavelength, flux vectors'
print,'* '
print,'* KEYWORDS'
print,'* ALLSPEC:'
print,'* FILE:'
print,'* SW, MW, LW:'
print,'* ICD:'
print,'* TITLE:'
print,' '
return
endif
if n_elements(root) eq 0 then root=''
if (not keyword_set(sw)) and (not keyword_set(mw)) and (not keyword_set(lw)) $
then allspec=1
if ifstring(file) then allspec=0
if not keyword_set(sw) then sw=0
if not keyword_set(mw) then mw=0
if not keyword_set(lw) then lw=0
if keyword_set(icd) then begin
; allspec=1
if not ifstring(icd) then icdfile='euve' else icdfile=icd
endif
if not keyword_set(title) then title=''
if keyword_set(allspec) then begin
sw=1 & mw=1 & lw=1
endif
nspec=fix(total((sw<1)+(mw<1)+(lw<1)))
if n_elements(file) gt 0 then nspec=1
;
for i=0,nspec-1 do begin
case 1 of
ifstring(file):
keyword_set(sw): begin
file='sw_spec'
sw=0
end
keyword_set(mw): begin
file='mw_spec'
mw=0
end
keyword_set(lw): begin
file='lw_spec'
lw=0
end
else: begin
print,' ERROR'
stop
return
end
endcase
print,' Reading ',file+'.fits'
dd=readfits(file+'.fits',h)
d=dd(*,0,0)
nax=getval('naxis',h)
case 1 of
nax eq 3: begin ;APEXTRACT
nax3=getval('naxis3',h)
t=getval('wat2_001',h)
k=strpos(t,'"')
t=strmid(t,k+1,80)
i1=0 & i2=0 & i3=0 & w0=0. & dw=0.
reads,t,i1,i2,i3,w0,dw
case 1 of
nax3 eq 4: begin ;variance weighting
sn=dd(*,0,3)
e=abs(d/sn)
end
nax3 eq 2: begin ;normal extraction
b=dd(*,0,1)
e=sqrt((b+d)*(b+d)+b*b)
end
else: begin
print,' NAXIS2=',nax3,' is invalid'
return
end
endcase
end
else: begin
w0=float(strtrim(getval(h,'crval1'),2))
dw=float(strtrim(getval(h,'cdelt1'),2))
e=dd(*,0,1)
end
endcase
nl=fix(strtrim(getval(h,'naxis1'),2))
lam=w0+findgen(nl)*dw
expt=float(getval(h,'exptime'))
if keyword_set(icd) then begin
head=intarr(400)
date=getval(h,'date',/noap)
head(10)=fix(strmid(date,3,2))
head(11)=fix(strmid(date,0,2))
head(12)=fix(strmid(date,6,2))
head(3)=80 ;EUVE type
head(5)=-fix(expt/60.) ;exp time in minutes
head(19)=10000
head(20)=fix(w0) & head(21)=fix(head(19)*(w0-fix(w0)))
head(22)=fix(dw) & head(23)=fix(head(19)*(dw-fix(dw)))
head(199)=333
k=strpos(file,']',0)
if k gt 0 then sf=strmid(file,k+1,32) else sf=file
kdat,icdfile,head,lam,d,e,-1,title=title+' '+sf,/islin
endif
endfor
;
if keyword_set(stp) then stop,'EUVESPEC>>>'
return
end
euvespec.pro;2 000640 000451 000017 00000000000 06615503115 017005 1euvespec.pro ustar 00fwalter users 000002 260633 ewlim.pro 000640 000451 000017 00000004603 06615503115 013601 0 ustar 00fwalter users 000002 260640 ;*****************************************************************************
PRO EWLIM,file,w1,w2,iprt=iprt,iint=iint
common comxy,xcur,ycur,zerr
;
if string(file) eq '-1' then begin
print,'*'
print,'* EWLIM - estimate minimun detectable equivalent widths'
print,'* calling sequence: EWLIM,file,w1,w2,iprt=iprt,iint=iint'
print,'* file : name of data file'
print,'* w1,w2: wavelength region for S/N determination (IINT not set)'
print,'* IPRT : set to print .LST file when done'
print,'* IINT : set for interactive processing'
return
endif
;
;
if n_params(0) eq 0 then read,' enter name of data file',file
if (n_params(0) lt 3) and (not keyword_set(iint)) then $
read,' enter wavelengths of continuum region',w1,w2
siglim=3.
;
if keyword_set(iint) then begin
read,' enter record number: ',irec
irec=fix(irec)
i0=irec & i1=irec
endif else begin
openw,lu,'ewlim.lst',/get_lun
printf,lu,' EWLIM output, run at ',systime(0)
printf,lu,' data file=',file
wlim=[w1<w2,w1>w2]
printf,lu,' wavelengths region = ',wlim
printf,lu,' '
i0=0 & i1=999
endelse
;
for i=i0,i1 do begin
print,i,i0,i1
gdat,file,h,w,f,e,i
if n_elements(h) lt 3 then goto,done
dw=w(1)-w(0)
title=strtrim(byte(h(100:160)>32),2)
f1=optfilt(f,e)
ff=f/f1
if keyword_set(iint) then begin
plot,w,ff
print,' use X command to reset limits, other to mark region'
blowup,-1
if (zerr eq 88) or (zerr eq 120) then begin
x1=xcur
blowup,-1
x2=xcur
!x.range=[x1<x2,x1>x2]
plot,w,ff
print,' Now mark region'
blowup,-1
endif
w1=xcur
blowup,-1
w2=xcur
wlim=[w1,w2]
endif
ii=fix(xindex(w,wlim)+0.5) ;tabinv,w,wlim,ii
np=1+ii(1)-ii(0)
tf=total(ff(ii(0):ii(1)))/float(np)
tk=(ff-tf)*(ff-tf)
tk=total(tk(ii(0):ii(1)))/float(np-1)
rms=sqrt(tk)
sn=abs(tf/rms)
minew=siglim/sn*3.*dw
minew=fix(minew*1000.+0.5)
si=string(i,format='(I3)')
if not keyword_set(iint) then $
printf,lu,'Record ',si,': min EW detectable (mA)=',minew,' Obj=',title
print,'Record ',si,': min EW detectable (mA)=',minew,' Obj=',title
endfor
done:
close,lu
free_lun,lu
if not keyword_set(iint) then begin
if keyword_set(iprt) then spawn_print,'ewlim.lst' else $
print,' output in EWLIM.LST'
endif
return
end
ewlim.pro;2 000640 000451 000017 00000000000 06615503115 015577 1ewlim.pro ustar 00fwalter users 000002 260640 ews.pro 000640 000451 000017 00000003623 06615503115 013266 0 ustar 00fwalter users 000002 260643 ;**********************************************************************
pro ews,file,r1,r2
ihlp=0
if n_params(0) eq 0 then file=-1
if not ifstring(file) then begin
if file eq -1 then ihlp=1
endif
if ihlp eq 1 then begin
print,' '
print,'* EWS - loop on EQWID through ICUR data file'
print,'* calling sequence: EWS,file,r1,r2'
print,'* file: name of .data file, or type for generic file'
print,'* r1,r2: records to measure, default=all'
print,'* output in EQWID.LST'
print,' '
return
endif
;
if not ifstring(file) then begin
if (file lt 0) or (file gt 110) then begin
print,' generic file type = ',file,' outside valid range'
return
endif
endif
if n_params(0) lt 3 then r2=9999
if n_params(0) lt 2 then r1=0
;
;set up eqwid params
;
print,' enter -1 for defaults'
read,' Enter wavelength of feature: ',lamcen
read,' Enter width for feature extraction: ',dl
read,' Enter wavelength of first background region: ',b1
read,' Enter wavelength of second background region: ',b2
read,' Enter width of first background region: ',db1
if db1 eq -1 then db2=-1 else read,' Enter width of second background region: ',db2
;
outfile='eqwid.lst'
openw,lu,outfile,/get_lun
printf,lu,' EQWID output, data file=',file,' run at ',systime(0)
printf,lu,' Feature at ',lamcen,' A, box width=',dl,' A.'
printf,lu,' Background bins begin at ',b1,b2,' A., box widths=',db1,db2
printf,lu,' '
printf,lu,'Record Equivalent width Object'
;
for i=r1,r2 do begin
gdat,file,h,w,f,e,i
if n_elements(h) eq 1 then goto,done ;record does not exist
ew=eqwid(w,f,lamcen,dl,b1,b2,db1,db2)
if ew gt -9998. then begin
sew=string(ew,format='(F6.2)')
printf,lu,i,sew,' ',strtrim(byte(h(100:158)),2)
print,i,sew,' ',strtrim(byte(h(100:158)),2)
endif ;-9999 if invalid region
endfor
;
done:
close,lu
free_lun,lu
print,' output in ',outfile
return
end
ews.pro;2 000640 000451 000017 00000000000 06615503115 014744 1ews.pro ustar 00fwalter users 000002 260643 export.txt 000640 000451 000017 00000022630 06615503116 014026 0 ustar 00fwalter users 000002 260650 To: IDL Users
From: Fred Walter
Date: 7/1/95
re: ICUR version 3.0 updates
See the end of the file for instructions on copying files.
***************************************************************************
June 1997
Upgraded to IDL 4.0.1. Removed all references to SET_VIEWPORT.
**********************************************************************
6/16/92
1. Fixed another bug with changing "frozen" parameters in the fitting
routine. This affected only frozen line widths in spectra with
non-linear wavelength scales.
2. Changed the operating philosophy. ICURSTARTUP.PRO now initializes the
common blocks and starts GOICUR. GOICUR has been modified to a procedure.
Read the documentation (updated 6/16/92). This change should be transparent
to most users. The only difference is to replace the .RUN GOICUR command
with GOICUR.
**********************************************************************
6/13/92
1. Fixed some nagging bugs in the fitting program. The errors are now more
realistic, and frozed parameters stay where they should.
2. When using data with real errors (epsilon code 30), the equivalent
width code now gives errors on the equivelent widths. COADD also adds
error bars now.
3. For GHRS users: option d does quick and dirty deconvolutions of
GHRS spectra. You will need access to the GHRS software, and will need
to put the PSF*.DAT files in ICURDATA.
- Thanks to R. Tarrall and J. Neff for uncovering more bugs hidden in dimly
lit corners of the code.
**********************************************************************
5/27/92
1. Finalized version 3.0
2. Revised documentation file ICUR.TEX
3. The .ICD files are now universally readable. Use ICD_VtoU to translate
files FTP'ed from VMS machines.
4. Cleared up some inconsistent variable names, and removed system variable
!debug.
**********************************************************************
5/24/92
1. Altered the way PLDATA knows how to plot bad data. If H(33)=30
(EPS vector is S/N vector), then the S/N vector is the bad data
vector, and can be plotted in PLDAT. Changed BDATA, PLDATA, and all
routines that call them. Note that this only affects PLDATA,3,... calls.
**********************************************************************
5/22/92
1. revised KDAT, GDAT, LDAT to check system architecture when .ICD
files are created. I can now read, under VMS, .ICD files written on a
Sparcstation. Further tests are needed to check the other direction.
2. Incorporated the real S/N vectors into the fitting
routines. Defined MODE=4 in ICFIT2.
3. altered LINEARWAVE to overwrite header word 199.
**********************************************************************
3/11/92
Added a spawn,'rm NL0:' at the end of ICUR for UNIX calls
**********************************************************************
3/6/92
1. QPLT: added switch for hard copies, /hc for print, hc=-1 to make .ps file
2. added /NOTITLE keyword to DAT_TO_ICD. Set it and you won't be asked
for titles.
3. common variable PSDEL, controlling the .PS output, has been revised.
- bit 0: set to save plot, 0 to delete
- bit 1: set to hold plot, 0 to plot
setting bit 1 overrides bit 0, i.e, psdel=2,3 are equivalent.
**********************************************************************
3/5/92
1. KDAT, GDAT, and LDAT replace KEEPDAT, GETDAT, LKDAT. Default extension
on all data files is .ICD. IDAT no longer exists as a data type
identifier. Standard files cannot be automatically identified, and
are prompted for at the first g command.
Use DAT_TO_ICD to translate data files. (GETDAT -> KDAT).
2. added w command to ICUR. This provides continuous readouts of cursor
position when it is inside the window. Procedure CWHERE.
3. Added USERDATA directory, at JEN's suggestion. Search path is:
- (1) Current directory, (2) USERDATA, (3) ICURDATA
USERDATA is read using GETENV, like ICURDATA
4. Removed KEEPADT and MKGENFIL from the library.
**********************************************************************
2/28/92
1. Fixed annoying bugs in FFCFIT: There was an itermittent 1/2 pixel error
in the line position. Also, the line width is now forced to be positive.
2. Revised logging option. Use the & command to turn it on and off within
ICUR. The log question in the main procedure has been dropped.
3. Simplified input to the main procedure. Just enter the file name.
No more -2 then the file name.
4. ICFIT: It now complains when you type an invalid option.
5. ICUR, FUN1: Now echo command to the X screen, and tell when its waiting
for more input.
Added procedure OPSTAT.
**********************************************************************
2/18/92
Defined a new data format to ensure compatability between systems. the
.ICD data files will replace the .DAT data files.
Added KDAT, GDAT, and LDAT to library. Added ancillary procedure LINEARWAVE.
**********************************************************************
2/7/92
Corrected factor of 2 error in Gausian Integrated Flux in FFLPFW.
**********************************************************************
1/31/92
Added new command - to permit freezing the Y scaling.
**********************************************************************
1/29/92
1. GOICUR has been modified to clean up the code.
2. An option to fix Y axis scaling has been added. Type _ (underscore)
to freeze Y scaling. The _ command toggles a flag. The flag is negated
by the Y or N commands.
3. The plotting options in ICFIT (option F) have been upgraded.
**************************************************************************
10/7/91 update
1. GOICUR.PRO, the main program, has been modified to permit use with
monochrome X-windows.
2. BBODY (option t) has been updated to permit use when the plot device
is set to the hard copy device (!). This lets the user overplot
black body curves on hard copies of spectra.
**************************************************************************
8/28/91 update
1. A bug in GHRSTOICUR was fixed. Spectra can again be cross correlated to
determine wavelength shifts.
2. GHRSTOICUR was modified to permit storage of summed direct downlink
data. Set the /dd keyword.
3. TABINV was updated using the GSFC version.
*************************************************************************
7/28/91
In response to complaints, the following improvements have now been incorporated
into ICUR.
1. The file UV.LIN should have been included for distribution. It is now.
The files HIGH.LIN and LOW.LIN are not used, and should be discarded.
2. All EXTRAC calls have been removed, as have some referrals to !sc3 and
!sc4 system variables. Further vestiges of version 1 will be removed as
necessary, and upon request.
3. A major problem with use of ICUR on TEK teminals has been corrected.
Basically, I forgot to put in a CASE statement for this plot device.
An improved cross correlation routine (a la Tonry and Davis), incorporating
estimates of rotational as well as radial velocities, will be added shortly.
Stay tuned.
****************************************************************************
****************************************************************************
****************************************************************************
How to copy ICUR to your home institution
All files are in 44156::rulupi$dka200:[users.public.icur]. You can access them
through the GUEST account. You can use VMS copy, i.e,
dir 44156"guest guest"::rulupi$dka200:[users.public.icur]*.* (to see what's there)
copy 44156"guest guest"::rulupi$dka200:[users.public.icur]xxx.xxx *
or FTP, i.e,
ftp sbast1.ess.sunysb.edu
cd rulupi$dka200:[users.public.icur]
ls (if you want to know what's there)
get ...
Note that the guest account is restricted - you cannot log into it to do
anything. You can do remote copies and directories, however.
The files you need are:
ICUR.COM : a sample command file for running ICUR.
ICUR.TLB : the text library containing all the procedures (1273 blocks).
ICUR.TEX : the LATEX documentation file (39+ pages).
ICURSTARTUP_V2.PRO: a main procedure which runs GOICUR.
ICUR.MSG : a file typed by GOICUR, containing messages.
*.LIN : 2 files (UV and OPT) containing incomplete line lists.
*.DAT : 5 files containing interstellar extinction curves. These are
needed for the u command. The files are:
NANDY.DAT, ORI.DAT, SAVMAT.DAT, SEATON.DAT, and SMC.DAT.
12 PSF*.DAT files, useful only for persons having GHRS
spectra and access to the GHRS software, for doing
quick-and-dirty deconvolutions.
SAMPLE.ICD: sample data file.
EXPORT.TXT: the list of current updates to the package (which you are now
reading).
optional files are:
FFIT2.FOR : the FORTRAN version of the fitting routine (optional).
FFIT2.EXE : The executable code for the FORTRAN fitting routine (optional).
read the documentation before you copy these files
Please report any problems to me at 44156::FWALTER or
FWALTER@SBAST1.ESS.SUNYSB.EDU. Also, please let me know if you actually use
the package, so that I can inform you about bug fixes.
*******************************************************************************
export.txt;2 000640 000451 000017 00000000000 06615503116 016247 1export.txt ustar 00fwalter users 000002 260650 factorial.pro 000640 000451 000017 00000001132 06615503116 014431 0 ustar 00fwalter users 000002 260637 ;**************************************************************************
function factorial,nn
n=nn
np=n_elements(nn)
k=where(nn gt 33,nk)
if nk gt 0 then begin
print,' FACTORIAL: WARNING: values >33 reset to 33'
n(nk)=33
endif
k=where(nn lt 1,nk)
if nk gt 0 then begin
n(nk)=1
endif
;
f=dblarr(np)+1.D0
for k=0,np-1 do begin
if n(k) gt 1 then for i=2L,long(n(k)) do f(k)=f(k)*i
endfor
if np eq 1 then begin
f=f(0)
case 1 of
f le 32767L: f=fix(f)
(f gt 32767L) and (f le 2147483647L): f=long(f)
else: f=float(f)
endcase
endif
return,f
end
factorial.pro;2 000640 000451 000017 00000000000 06615503116 017264 1factorial.pro ustar 00fwalter users 000002 260637 fchisq.pro 000640 000451 000017 00000001507 06615503116 013751 0 ustar 00fwalter users 000002 260647 ;*******************************************************************
FUNCTION FCHISQ,nfree,y,yfit,weight
; called by FFIT2
; returns chi-2 value
;
if n_params(0) lt 3 then begin
print,' '
print,' * FCHISQ'
print,' * returns chi-2 of fit YFIT to data Y, given errors SIG'
print,' * calling sequence: CHI2=FCHISQ(NFREE,Y,YFIT,WEIGHT)'
print,' * NFREE: number of free parameters (greater than 0)'
print,' * Y: data'
print,' * YFIT: fit to data'
print,' * WEIGHT: weights of data, default = uniform weighting'
print,' '
return,-99. ;not enough parameters
endif
;
IF NFREE le 0 then return,-1. ; not enough degrees of freedom
if n_params(0) lt 4 then weight=y*0.+1. ;equal weights
CHISQ=WEIGHT*(Y-YFIT)*(Y-yfit)
FCHI=total(CHISQ)/float(nFREE)
RETURN,fchi
END
fchisq.pro;2 000640 000451 000017 00000000000 06615503116 016107 1fchisq.pro ustar 00fwalter users 000002 260647 ffcfit.pro 000640 000451 000017 00000004607 06615503117 013735 0 ustar 00fwalter users 000002 260660 ;************************************************************************
pro ffCFIT,NTR,Mode,chisq,iflag,cii=cii
; called by ffit2
;C** CFIT - FORMERLY FWCFIT
;C** CALLED BY FFIT; CALLS FFCURFIT
;C** INPUTS: N=NUMBER OF POINTS
;C** NTR = NUMBER OF FREE PARAMETERS
; version of ffcfit to call fortran version of FFCURFIT
COMMON CFT,X0,Y0,SIG0,epsvect
COMMON CURVE,A,EA,IFIXT
COMMON PLCR,YF
common fortflag,fflg
common custompars,dw,lp,x2,x3,x4,x5
;
if n_elements(fflg) eq 0 then fflg=0
chisq=0.
n=where(ifixt gt 0)
if n(0) ne -1 then nvar=n_elements(n) else nvar=0
if nvar le 0 then begin
print,' CFIT error: Number of variables =',nvar,' - returning to main'
return
endif
;
n=n_elements(x0)
nfree=n-nvar ;degrees of freedom
if nfree le 0 then return
;
x=double(x0)
y=double(y0)
sig=double(sig0)
a=double(a)
np=fix(n_elements(a))
nlines=(np-3)/3
for i=1,nlines do a(i*3+1)=a(i*3+1)-0.5 ;calibration of binning offset
NTERMS=fix(NTR)
IF NTerms LT 3 then NTerms=3
ntr=nterms
ishape=fix(nlines+1)
FLAMDA=.001D0
CHIOLD=0.D0
nt=nterms
case 1 of
mode lt 0: weight=1./abs(y) ;1/y weighting
mode eq 0: weight=fltarr(n)+1. ;uniform weighting
else: begin
weight=1./sig/sig ;weighting by variance
kbad=where(sig le -9998,nbad)
if nbad gt 0 then weight(kbad)=0.
end
endcase
if fflg gt 0 then begin ;fixed array sized for FORTRAN calls
x1=fltarr(750) & x1(0)=float(x0)
y1=fltarr(750) & y1(0)=float(y0)
w1=fltarr(750) & w1(0)=float(weight)
a1=fltarr(18) & a1(0)=a
ea1=fltarr(18) ; & ea1(0)=ea
ifixt1=intarr(18) & ifixt1(0)=ifixt
flamda=.001
chisq1=0.
iflag=0
endif
;
;*** ;special setups
ni=n_elements(ifixt)
if keyword_set(cii) then begin
if ni ge 8 then ifixt(7)=0
if ni ge 14 then ifixt(13)=0
endif
;***
;
for IGLOO=0,12 do begin
case 1 of
fflg eq 2: calltest,'FCURFIT2','FCURFIT2_EXE', $
x1,y1,w1,a1,ea1,ishape,ifixt1,np,nt,flamda,chisq1,iflag ;fortran version
fflg eq 1: fCURFIT,x1,y1,w1,a1,ea1,ishape,ifixt1,np,nt,flamda,chisq1,iflag ;fortran version
else: FFCURFIT,Nfree,weight,FLAMDA,CHISQ,IFLAG,cii=cii ;IDL version
endcase
IF IFLAG NE 0 then return
IF(ABS(CHIOLD-CHISQ) LT .0003)then GOTO,RET
CHIOLD=CHISQ
endfor
;
ret:
a=float(a)
ea=float(ea)
for i=1,nlines do a(i*3+1)=a(i*3+1)+0.5 ;calibration of binning offset
return
end
ffcfit.pro;2 000640 000451 000017 00000000000 06615503117 016053 1ffcfit.pro ustar 00fwalter users 000002 260660 ffcurfit.pro 000640 000451 000017 00000004756 06615503117 014314 0 ustar 00fwalter users 000002 260663 ;**********************************************************************
pro FFCURFIT,Nfree,weight,FLAMDA,CHISQ1,IFLAG,cii=cii
; called by ffcfit in ffit2
;C** CURFIT FROM BEVINGTON VIA FWCFIT
;C** TYPED IN 6/25/82
;C** modified 7/2/84 to permit freezing of parameters, following
;C** M. Lampton's WHIZ algorithm
COMMON CFT,X,Y,SIGMAY,eeeeee
COMMON PLCR,YFIT
COMMON CURVE,A,EA,IFIXT
common custompars,dw,lp,x2,x3,x4,x5
common ffits,lu4
if n_elements(lu4) eq 0 then lu4=-1
ea=a*0.
iflag=-1
TW=TOTAL(WEIGHT)
;
n=where(ifixt gt 0,NVAR)
INDX=FIX(N)
;
BETA=dblarr(nvar,/nozero)
ALPHA=dblarr(nvar,nvar,/nozero)
array=dblarr(nvar,nvar,/nozero)
FUN=FUNGUS(x,A,cii=cii)
deriv=GDERIV(x,a)
mm=min(abs(deriv)>0.)
isign=deriv/(abs(deriv)>mm)
logd=alog10(abs(deriv))
logw=alog10(weight+1.E-9)
z=weight*(y-fun)
for j=0,nvar-1 do beta(j)=total(z*deriv(*,indx(j)))
;
for j=0,nvar-1 do begin
JI=INDX(J)
for k=0,j DO begin
KI=INDX(K)
z=(logW+logD(*,JI)+logD(*,KI))>(-37.)
ALPHA(J,K)=TOTAL(isign(*,ji)*isign(*,KI)*10^(z))
endfor
endfor
;
yfit=fungus(x,a,cii=cii)
IFLAG=-1
for J=0,NVAR-1 do begin
for K=0,J do ALPHA(K,J)=ALPHA(J,K)
IF ALPHA(J,J) EQ 0. then goto,curerr
endfor
IFLAG=0
CHISQ1=FCHISQ(NFREE,y,yfit,weight)
lab71:
sqralpha=sqrt(ABS(alpha))
for J=0,NVAR-1 do begin
dj=sqrt(alpha(j,j))
for K=0,NVAR-1 do ARRAY(J,K)=ALPHA(J,K)/dj/sqrALPHA(K,K)
ARRAY(J,J)=1.D0+FLAMDA
endfor
;
iflag=-99
arr=array
array=invert(arr)
;
b=a
dal=dblarr(nvar,/nozero)
for j=0,nvar-1 do dal(j)=alpha(j,j)
dal=sqrt(dal)
for J=0,NVAR-1 do begin
JI=INDX(J)
b(ji)=b(ji)+total(beta*array(j,*)/dal(j)/dal)
;for K=0,NVAR-1 do B(JI)=b(ji)+BETA(K)*ARRAY(J,K)/SQRT(ALPHA(J,J))/sqrt(ALPHA(K,K))
endfor
YFIT=FUNGUS(x,B,cii=cii)
CHISQR=FCHISQ(NFREE,y,yfit,weight)
case 1 of
chisq1 lt chisqr: begin
FLAMDA=10.D0*FLAMDA
IF FLAMDA LT 101.D0 then GOTO,lab71
IFLAG=-92
RETURN
end
chisq1 ge chisqr: begin
a=b
saa=SQRT(ABS(ARRAY/ALPHA))
for J=0,NVAR-1 do begin
JI=INDX(J)
; EA(JI)=SQRT(ARRAY(J,J)/ALPHA(J,J))
EA(JI)=SAA(J,J)
endfor
end
endcase
FLAMDA=FLAMDA/10.D0
CHISQ1=CHISQR
iflag=0
RETURN
;
curerr:
bell
print,' *** CURFIT failed to converge'
if lu4 le 0 then return
printf,lu4,'* CURFIT ERRORS: IFLAG,J,Nvar=',IFLAG,J,Nvar
printf,lu4,'* I,A(I),DERIV(I),BETA(I),ALPHA(I,J=(1,N))'
for I=0,Nvar-1 do begin
ji=indx(i)
printf,lu4,'* ',I,A(jI),DERIV(jI),BETA(I),(ALPHA(I,0:nvar-1))
endfor
RETURN
END
ffcurfit.pro;2 000640 000451 000017 00000000000 06615503117 016774 1ffcurfit.pro ustar 00fwalter users 000002 260663 fffitstat.pro 000640 000451 000017 00000006073 06615503117 014466 0 ustar 00fwalter users 000002 260670 ;***************************************************************************
pro fffitstat,chisq,iflag,bsca
COMMON CFT,X,Y,SIG,E
COMMON CURVE,A,EA,IFIXT
COMMON PLCR,YF
common ffits,lu4
common icurunits,xunits,yunits,title
if n_elements(lu4) eq 0 then lu4=-1
n=n_elements(x)
nterms=n_elements(a)
nlines=(nterms-3)/3
L1=1
CASE 1 OF
nterms le 3: printf,lu4,'* No lines fit'
(NTERMS gt 3) and (nterms le 6): begin ;1 line
printf,lu4,'* BEST FIT LINE INTENSITY=', $
string(a(3),'(G10.3)'),' +/-',string(ea(3),'(g10.3)')
print,' BEST FIT LINE INTENSITY=',a(3),' +/-',Ea(3)
IF NTERMS ge 5 then begin
printf,lu4,'* FIT BIN =',a(4),' +/-',ea(4),' ; GIVEN BIN',bsca
print,' FIT BIN =',a(4),' +/-',ea(4),' ; GIVEN BIN',bsca
IF NTERMS ne 5 then begin
printf,lu4,'* FIT HWHM=',a(5),' +/-',ea(5),' BINS'
print,' FIT HWHM=',a(5),' +/-',ea(5),' BINS'
endif ;ne 5
endif ;ge 5
end
nterms gt 6: begin ;2 or more lines
while nlines ge 2 do begin
l2=l1+1
k1=l1*3
k2=l2*3
printf,lu4,'* LINES',l1,' AND',l2
printf,lu4,'* LINE AMPLITUDES:', $
string(a(k1),'(G10.3)'),' +/-',string(ea(k1),'(g10.3)'),' ', $
string(a(k2),'(G10.3)'),' +/-',string(ea(k2),'(g10.3)')
printf,lu4,'* LINE POSITIONS :', $
string(a(k1+1),'(G10.3)'),' +/-',string(ea(k1+1),'(g10.3)'),' ', $
string(a(k2+1),'(G10.3)'),' +/-',string(ea(k2+1),'(g10.3)')
printf,lu4,'* LINE HWHMS :', $
string(a(k1+2),'(G10.3)'),' +/-',string(ea(k1+2),'(g10.3)'),' ', $
string(a(k2+2),'(G10.3)'),' +/-',string(ea(k2+2),'(g10.3)')
nlines=nlines-2
l1=l1+2
endwhiLE
if nlines eq 1 then begin
l2=l1+1
k1=l1*3
k2=l2*3
printf,lu4,'* LINE ',l1
printf,lu4,'* LINE AMPLITUDE:',a(k1),' +/-',ea(k1)
printf,lu4,'* LINE POSITION:',a(k1+1),' +/-',ea(k1+1)
printf,lu4,'* LINE HWHM:',a(k1+2),' +/-',ea(k1+2)
endif
end
else:
endcase
;
schs=string(chisq,'(F8.3)')
printf,lu4,'* CHI-SQUARE per degree of freedom:',schs
CHISQ=CHISQ*FLOAT(N-NTERMS)
sflag=string(iflag,'(I3)')
sn=string(n,'(I4)')
snt=string(nterms,'(I3)')
schs=string(chisq,'(F8.3)')
printf,lu4,'* Flag= ',sflag,' Chi**2=',schs,' Data points= ',sn,' Free parameters= ',snt
print,'* Flag= ',sflag,' Chi**2=',schs,' Data points= ',sn,' Free parameters= ',snt
;
NPL=6
imax=fix((nterms-1)/NPL)
SNPL=STRTRIM(NPL,2)
SF='G11.4'
fmt1='("",A,"A:",'+SNPL+'(" ",'+SF+'))'
fmt2='(" EA:",'+SNPL+'(" ",'+SF+'))'
for i=0,imax do begin
s=string(i*npl,'(i2)')+'-'+string((i*npl+npl-1)<(nterms-1),'(i2)')+' '
printf,lu4,S,a(i*NPL:(nterms-1)<(i*NPL+NPL-1)),format=fmt1
printf,lu4,ea(i*NPL:(nterms-1)<(i*NPL+NPL-1)),format=fmt2
; printf,lu4,'* ',S,' A:', a(i*NPL:(nterms-1)<(i*NPL+NPL-1))
; printf,lu4,'* EA:',ea(i*NPL:(nterms-1)<(i*NPL+NPL-1))
endfor
;
RETURN
END
fffitstat.pro;2 000640 000451 000017 00000000000 06615503117 017336 1fffitstat.pro ustar 00fwalter users 000002 260670 fffixlines.pro 000640 000451 000017 00000001430 06615503117 014626 0 ustar 00fwalter users 000002 260666 ;***********************************************************************
pro fffixlines,ncam
COMMON CFT,X,Y,SIG,E
COMMON CURVE,A,EA,IFIXT
common ffits,lu4
np=n_elements(a)
nlines=(np-3)/3
if nlines le 0 then return
for i=0,nlines-1 do begin ;are lines to be fixed here?
klin=(i+1)*3+1
if a(klin) lt 0 then begin
IF NCAM EQ 10 then A(klin+1)=3.7 ;IIDS
IF (NCAM EQ 20) OR (NCAM ge 30) then A(klin+1)=2.0
A(klin)=ABS(A(klin))
wl=lint(x,A(klin))
print,' LINE AT ',wl,' A. ENTER proper wavelength, 0 to go on'
READ,wnew
if wnew gt 0. then begin
Da=(WNEW-WL)/apb
A(klin)=A(klin)+DA
wl=lint(x,a(klin))
print,' LINE FIXED AT BIN',a(klin),' WAVELENGTH=',wl
endif
endif
endfor
return
end
fffixlines.pro;2 000640 000451 000017 00000000000 06615503117 017651 1fffixlines.pro ustar 00fwalter users 000002 260666 ffile.pro 000640 000451 000017 00000001167 06615503120 013560 0 ustar 00fwalter users 000002 260667 ;*************************************************************************
function ffile,disk,direc,mfile,ext,notify=notify
if n_params(0) lt 1 then return,0
if n_elements(disk) eq 0 then return,0
case 1 of
n_params(0) eq 1: file=disk
n_params(0) eq 2: file=disk+'.'+direc
n_params(0) eq 3: file=disk+direc+'.'+mfile
else: file=disk+':'+direc+mfile+'.'+ext
endcase
on_ioerror,dne
get_lun,lu
openr,lu,file
igo=1 ;file exists
close,lu
goto,ret
dne: igo=0
if keyword_set(notify) then begin
bell
print,' *** FILE ',file,' was not found. ***'
endif
on_ioerror,null
ret: free_lun,lu
return,igo
end
ffile.pro;2 000640 000451 000017 00000000000 06615503120 015524 1ffile.pro ustar 00fwalter users 000002 260667 ffinflux.pro 000640 000451 000017 00000002533 06615503120 014310 0 ustar 00fwalter users 000002 260674 ;****************************************************************************
pro FFINFLUX,igo
; called by FFIT2
; determines integrated fluxes
COMMON CFT,X,Y,SIG,E
common ffits,lu4
if n_elements(lu4) eq 0 then lu4=-1
imax=n_elements(x)
if n_params(0) eq 0 then igo=-1 ;igo not equal to -1 to query for fluxes
I1=0
I2=IMAX-1
WD=(X(IMAX-1)-X(0))/FLOAT(IMAX-1)
;
restart:
k=indgen(i2-i1+1)+i1
wsat=where(e lt -1750)
nsat=n_elements(wsat)
if wsat(0) eq -1 then nsat=0
flux=total(y(k))
sk=sig(k)
kbad=where(sk ge 998.)
if kbad(0) ne -1 then sk(KBAD)=0.
fer=total(sk*sk)
W1=X(I1)
W2=X(I2)
DW=W2-W1+WD
FLUX=FLUX*DW/FLOAT(I2-I1+1)
FER=SQRT(FER)*DW/FLOAT(I2-I1+1)
sw1=string(w1,'(F9.3)')
sw2=string(w2,'(F9.3)')
if lu4 gt 0 then printf,lu4,'* Flux from',sW1,' to',SW2,' Ang. =',flux,' +/-',fer
print,' Flux from',SW1,' to',SW2,' Ang. =',flux,' +/-',fer
IF NSAT gt 0 then begin
if lu4 gt 0 then printf,lu4,nsat,'* SATURATED PIXELS INCLUDED'
print,nsat,' SATURATED PIXELS INCLUDED'
endif
;
if igo eq -1 then return
;
print,' WAVELENGTH RANGE ACCESSIBLE IS ',x(0),' TO ',x(imax-1)
print,' ENTER RANGE OF INTEREST'
READ,W1
IF W1 LE 0. then RETURN
READ,W2
IF W2 LE W1 then begin
t=w2
w2=w1
w1=t
endif
W1=w1>X(0)
W2=w2<X(IMAX-1)
i1=fix(xindex(x,w1)+0.5) ;tabinv,x,w1,i1
i2=fix(xindex(x,w2)+0.5) ;tabinv,x,w2,i2
goto,restart
return
END
ffinflux.pro;2 000640 000451 000017 00000000000 06615503120 017012 1ffinflux.pro ustar 00fwalter users 000002 260674 ffit.fit 000640 000451 000017 00000004223 06615503120 013365 0 ustar 00fwalter users 000002 260700
**************************************************************
* FFIT, version 8.0, run on Tue Sep 16 1997 at 12:34
*
* IMAGE: 0
* Background bins used: 50 - 445
* BACKGROUND FIT: BEGINNING ESTIMATES: -24.7321 0.000000 0.000000
* LINE # 1: BEGINNING ESTIMATES: 40.5134 221.819 5.00000
* LINE # 2: BEGINNING ESTIMATES: 44.5313 230.703 5.00000
* Reciprocal Dispersion: 2.6383 Angstroms per bin
* LINES 1 AND 2
* LINE AMPLITUDES: -26.7 +/- 2.53 62.1 +/- 1.09
* LINE POSITIONS : 170. +/- 0.500 198. +/- 0.605
* LINE HWHMS : 10.9 +/- 1.16 43.5 +/- 1.14
* CHI-SQUARE per degree of freedom:4596.092
* Flag= 0 Chi**2=******** Data points= 396 Free parameters= 9
0- 5 A: 1.436 -0.03466 8.655E-05 -26.68 170.1 10.86
EA: 0.4689 0.003132 8.596E-06 2.529 0.4998 1.160
6- 8 A: 62.11 198.4 43.46
EA: 1.093 0.6051 1.137
* Instrumental width= 5.00 bins= 13.225 Angstroms
* Line 1 at 1104.945+/- 1.318 Angstroms
* Gaussian sigma of line is 28.734 Angstroms; FWHM=33.820 Angstroms
* Deconv. width = 31.127+/- 6.277 Angstroms; equiv. Vsin i=8445.3km/s
* Flux in line 1 + background: 716. +/- 4.81
* Extrapolated background flux: 1.68e+03 +/- 1.24
* Net Flux, line 1 : -967. +/- 4.97 erg/s/cm2
* Integrated Gaussian Flux= -960.928 erg/s/cm2; G/net= 0.99
* EW of line 1 is 48.469 +/- 0.003 Angstroms
*--
* Line 2 at 1179.435+/- 1.597 Angstroms
* Gaussian sigma of line is ****** Angstroms; FWHM=****** Angstroms
* Deconv. width =134.633+/- 6.156 Angstroms; equiv. Vsin i=******km/s
* Flux in line 2 + background: 7.07e+03 +/- 16.3
* Extrapolated background flux: -1.43e+03 +/- 1.24
* Net Flux, line 2 : 8.50e+03 +/- 16.3 erg/s/cm2
* Integrated Gaussian Flux= 8947.24 erg/s/cm2; G/net= 1.05
* EW of line 2 is 3373.353 +/- 0.012 Angstroms
*--
* Flux from 656.071 to 1698.185 Ang. = 7022.02 +/- 17.0332
ffit2.exe 000640 000451 000017 00000045000 06615503120 013445 0 ustar 00fwalter users 000002 260710 ¨ 0 D X 0205