Commit 68eac64f authored by Monica Rainer's avatar Monica Rainer
Browse files

Delete straight_giano_2D_v1_1_old.f

parent 488f26b1
program straight_giano_2D
C
C Compile : f77 -o straight_giano_2D {name_of_this_file} -lcfitsio
C
C Synopsis: straight_giano_2D in_file_name out_file_name [options]
C
C Purpose: eliminate curvature of orders and slit-tilt in echellogram
C
C Input (from line command, see above):
C - The name of a R*4 fits file with the echellogram to be Y-straightned
C
C - Options:
C I=0 default: does not interpolate i.e. shifts each pixel
C belonging to an order by an integer number of pixels
C I=1 linearly interpolate pixels in X and Y
C
C T=1 default: correct for slit tilt
C T=0 does not correct for slit tilt
C
C DY=nn Correction for rigid shift along Y (default=0)
C
C DT=tt Correction for rigid tilt of slit (default=0)
C
C V=xx Value of pixels at order separation (every 41 pixels)
C default=0
C
C Output
C - The output fits file
C
implicit integer(i-n)
implicit real (a-h,o-z)
dimension rima_in(2048,2048),rima_out(2048,2048)
dimension naxes(2)
character filein*200,fileout*200
character card*400
logical simple,extend,anyf,interp
C coefficients defining slit-tilt (linear with orders)
dimension slit_tilt(50)
C coefficients defining Y-deformation for each order,
C 50 orders from 32 to 81
dimension Xmin(50),C0(50),C1(50),C2(50)
C
data C0/02,68,130,188,242,294,343,390,435,478,520,561,600,639,
!!! orders 32 33 34 35 36 37 38 39 40 41 42 43 44 45
, 676,713,750,786,821,857,891,926,960,994,1029,1063,
!!! orders 46 47 48 49 50 51 52 53 54 55 56 57
, 1097,1132,1167,1201,1236,1271,1307,1343,1378,1414,
!!! orders 58 59 60 61 62 63 64 65 66 67
, 1451,1488,1526,1564,1603,1642,1681,1721,1762,1803,
!!! orders 68 69 70 71 72 73 74 75 76 77
, 1845,1887,1930,1974/
!!! orders 78 79 80 81
C iDYorder= width of each order (pixels)
data iDYorder/40/
C
C analytical expressions for deformation parameters
do iorder=32,81
io=iorder-32+1
Xord=float(iorder)
Xmin(io)=833.2-24.51*Xord+230.3*sqrt(Xord-25.0)
C1(io)=0
C2(io)=120.5+0.765*Xord
enddo
C fixed parameters for fitsio routines
data istatus,iblocksize,naxis/0,1,2/
data nullval,igroup,ifpix,maxdim/0,1,1,2/
data simple,extend,anyf/.true.,.true.,.false./
C useful parameters
data lunit_in/21/ ! unit to read fits files
data lunit_out/22/! unit to write output fits file
call getarg(1,filein)
call getarg(2,fileout)
if(filein(1:1).eq." " .or. fileout(1:1).eq." ") then
write(6,'(a)')
+ 'Usage: sraight_giano_2D in_file_name out_file_name [options]'
write(6,'(a)')
+ 'Options are: '
write(6,'(a,a)')
+ 'I=0 default: no interpolation i.e. shifts each pixel',
+ ' belonging to an order by an integer number of pixels.'
write(6,'(a)')
+ 'I=1 linearly interpolate values of pixels'
write(6,'(a,a)')
+ 'T=1 default: corrects for slit-tilt'
write(6,'(a)')
+ 'T=0 does not correct for slit-tilt'
write(6,'(a)')
+ 'DY=nn Correction for rigid shift along Y (default=0)'
write(6,'(a)')
+ 'DT=tt Correction for rigid tilt of slit (default=0)'
write(6,'(a,a)')
+ 'V=xx Value of pixels at order separation',
+ ' (every 41 pixels), default=0'
stop
endif
C Set default values
Val_separator=0
Tilt_apply=1.0
Shift_Y=0.0
Shift_tilt=0.0
interp=.false.
C Read options
do iarg=3,10
call getarg(iarg,card)
if(card(1:1).eq." ") goto 10
if(card(1:3).eq."I=1") interp=.true.
if(card(1:3).eq."I=0") interp=.false.
if(card(1:3).eq."T=1") Tilt_apply=1.0
if(card(1:3).eq."T=0") Tilt_apply=0.0
if(card(1:3).eq."DY=") read(card(4:80),*,err=995) Shift_Y
if(card(1:3).eq."DT=") read(card(4:80),*,err=995) Shift_tilt
if(card(1:2).eq."V=") read(card(3:80),*,err=995) Val_separator
enddo
10 continue
C
C analytical expression for slit-tilt: varies linearly with Ypos of order
data ST0,ST1/-0.0240,-0.137/
do iorder=32,81
io=iorder-32+1
Yord_norm=(C0(io)-1024.)/2048.
slit_tilt(io)=ST0+ST1*Yord_norm
slit_tilt(io)=slit_tilt(io)+Shift_tilt
slit_tilt(io)=slit_tilt(io)*Tilt_apply
enddo
istatus=0
irwmode=0 ! readonly
call ftopen(lunit_in,filein,irwmode,iblocksize,istatus)
if(istatus.ne.0) go to 994
C Get primary header parameters
call ftghpr(lunit_in,maxdim,simple,
, ibitpix,naxis,naxes,ipcount,igcount,extend,istatus)
C copy matrix into R*4 vector
nelem=naxes(1)*naxes(2)
call ftgpve(lunit_in,igroup,ifpix,nelem,nullval,
, rima_in,anyf,istatus)
C Create output image
ibitpix=-32 ! Real numbers
irwmode=1 ! =1 for readwrite
! open scratch file and rename it afterwards
call system('rm -f testscr_Y_out')
call ftinit(lunit_out,'testscr_Y_out',iblocksize,istatus)
call ftphpr
+ (lunit_out,simple,ibitpix,naxis,naxes,0,1,extend,istatus)
C copy keywords to output image
call copy_hdr_key(lunit_in,lunit_out)
C add keywords for straight_2D reduction
call ftpkys(lunit_out,"STR2D_P1",filein(1:60),
, 'straight_2D param',istatus)
call ftpkys(lunit_out,"STR2D_P2",fileout(1:60),
, 'straight_2D param',istatus)
if(interp) then
write(card(1:3),'(a)') "I=1"
else
write(card(1:3),'(a)') "I=0"
endif
call ftpkys(lunit_out,"STR2D_P3",card(1:3),
, 'straight_2D param',istatus)
write(card(1:13),'(a,f10.2)') "DY=",Shift_Y
call ftpkys(lunit_out,"STR2D_P4",card(1:13),
, 'straight_2D param',istatus)
write(card(1:13),'(a,f10.2)') "DT=",Shift_Tilt
call ftpkys(lunit_out,"STR2D_P5",card(1:13),
, 'straight_2D param',istatus)
write(card(1:13),'(a,1pe10.2)') "V=",Val_separator
call ftpkys(lunit_out,"STR2D_P6",card(1:13),
, 'straight_2D param',istatus)
C close input image
call ftclos(lunit_in,istatus)
C initialize output image
do iy=1,2048
do ix=1,2048
rima_out(ix,iy)=0.0
if(mod(iy,41).eq.0)
+ rima_out(ix,iy)=Val_separator ! to show order edges
enddo
enddo
C straighten order-by order
do iorder=32,81
io=iorder-32+1
jYout=1+(iorder-32)*41
do iX=1,2048
XN=(float(iX)-Xmin(io))/2048.
Ystart=C0(io)+C1(io)*XN+C2(io)*XN*XN
Ystart=Ystart+Shift_Y+2 ! Value corrected 17/11/2016
! tricks to better Y-centring lowest orders
if(iorder.lt.63) Ystart=Ystart-1
if(iorder.lt.36) Ystart=Ystart-1
do jY=1,iDYorder
Yval=Ystart+float(jY-1)-1
dY=float(jY)-20.
Xval=float(iX)+dY*slit_tilt(io)
kY=jYout+jY-1
if(kY.ge.1 .and. kY.le.2048) rima_out(iX,kY)=
= finterp(interp,Xval,Yval,rima_in,2048)
enddo
enddo
enddo
C
C Write processed image in output fits file
C
222 continue
nelem=naxes(1)*naxes(2)
call ftppre(lunit_out,1,1,nelem,rima_out,istatus)
call ftclos(lunit_out,istatus)
write(card,'(2a)') 'mv -f testscr_Y_out ',fileout
call system(card)
go to 999
C
994 continue
write(6,*) 'file',filein,' cannot be opened'
go to 999
995 continue
write(6,*) "Wrong option"
write(6,'(a)') card
go to 999
999 continue
! remove output scratch file if previously created
call system('rm -f testscr*')
end
real function finterp(interp,Xval,Yval,rima_in,NXY)
implicit integer(i-n)
implicit real (a-h,o-z)
logical interp
dimension rima_in(NXY,NXY)
C Avoid extrapolation
if(Yval.gt.2048.5 .or. Yval.lt.-0.5 .or.
+ Xval.gt.2048.5 .or. Xval.lt.-0.5) then
finterp=0.
return
endif
iY1=nint(Yval-0.5)
iY2=nint(Yval+0.5)
iY1=max(iY1,1)
iY2=max(iY2,1)
iY1=min(iY1,NXY)
iY2=min(iY2,NXY)
if(interp) then
dY=Yval-float(IY1)
else
dY=0.0
endif
iX1=nint(Xval-0.5)
iX2=nint(Xval+0.5)
iX1=max(iX1,1)
iX2=max(iX2,1)
iX1=min(iX1,NXY)
iX2=min(iX2,NXY)
if(interp) then
dX=Xval-float(IX1)
else
dX=0.0
endif
finterp=rima_in(IX1,IY1)*(1-dY)*(1-dX)+
+ rima_in(IX1,IY2)*dY*(1-dX)+
+ rima_in(IX2,IY1)*(1-dY)*dX+
+ rima_in(IX2,IY2)*dY*dX
return
end
subroutine copy_hdr_key(lunit_in,lunit_out)
!
! Copy and creates keywords in headers in output frame
!
! nima : # of image in input file list, used to decide if/what to copy
! from input frame
!
! lunit_in : logical unit number of input fits file
! lunit_out : logical unit number of output fits file
!
implicit integer(i-n)
implicit real(a-h,o-z)
character char*80
istatus=0
call ftghsp(lunit_in,nkeys,nspace,istatus) ! get number of keywords
do j=1,nkeys
call ftgrec(lunit_in,j,char,istatus) ! read j-th record keyword
if(char(1:8).eq.'SIMPLE ') goto 10 ! do not copy this keyword
if(char(1:8).eq.'BITPIX ') goto 10
if(char(1:8).eq.'NAXIS ') goto 10
if(char(1:8).eq.'NAXIS1 ') goto 10
if(char(1:8).eq.'NAXIS2 ') goto 10
if(char(1:8).eq.'EXTEND ') goto 10
if(char(1:8).eq.'BZERO ') goto 10
if(char(1:8).eq.'BSCALE ') goto 10
!if(char(1:8).eq.'COMMENT ') goto 10
!if(char(1:8).eq.'IMANAME ') goto 10
!if(char(1:8).eq.'IMAALIAS') goto 10
!if(char(1:8).eq.'GROUP_N ') goto 10
!if(char(1:8).eq.'GROUP_I ') goto 10
!if(char(1:8).eq.'MULTND_I') goto 10
!if(char(1:8).eq.'MULTI_OP') goto 10
!if(char(1:8).eq.'MULTICMD') goto 10
!if(char(1:8).eq.'TIME ') goto 10
!if(char(1:8).eq.'DIT ') goto 10
!if(char(1:8).eq.'EXPTIME ') goto 10
!if(char(1:8).eq.'PROGTINT') goto 10
call ftprec(lunit_out,char,istatus) ! write j-th record keyword
10 continue
enddo
return
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment