straight_giano_2D_v1_1_old.f 10.1 KB
Newer Older
Monica Rainer's avatar
Monica Rainer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
	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