1*======================================================================= 2* 3* WCSLIB 7.7 - an implementation of the FITS WCS standard. 4* Copyright (C) 1995-2021, Mark Calabretta 5* 6* This file is part of WCSLIB. 7* 8* WCSLIB is free software: you can redistribute it and/or modify it 9* under the terms of the GNU Lesser General Public License as published 10* by the Free Software Foundation, either version 3 of the License, or 11* (at your option) any later version. 12* 13* WCSLIB is distributed in the hope that it will be useful, but WITHOUT 14* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 15* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16* License for more details. 17* 18* You should have received a copy of the GNU Lesser General Public 19* License along with WCSLIB. If not, see http://www.gnu.org/licenses. 20* 21* Author: Mark Calabretta, Australia Telescope National Facility, CSIRO. 22* http://www.atnf.csiro.au/people/Mark.Calabretta 23* $Id: tspc.f,v 7.7 2021/07/12 06:36:49 mcalabre Exp $ 24*======================================================================= 25 26 PROGRAM TSPC 27*----------------------------------------------------------------------- 28* 29* TSPC tests the spectral transformation driver routines for closure. 30* 31*----------------------------------------------------------------------- 32* Maximum length of spectral axis - see CLOSURE. 33 INTEGER NSPEC 34 PARAMETER (NSPEC = 8001) 35 36 INTEGER CLOSURE, NAXISJ, NFAIL 37 DOUBLE PRECISION C, CDELTX, CRPIXJ, CRVALX, MARS(0:6), RESTFRQ, 38 : RESTWAV, X1, X2 39 40 COMMON /SPECTRO/ MARS 41 42 DATA C /2.99792458D8/ 43 44* KPNO MARS spectrograph grism parameters. 45 DATA MARS /4.5D5, 1D0, 27D0, 1.765D0, -1.077D6, 3D0, 5D0/ 46*----------------------------------------------------------------------- 47 WRITE (*, 10) 48 10 FORMAT ('Testing closure of WCSLIB spectral transformation ', 49 : 'routines (tspc.f)',/, 50 : '--------------------------------------------------', 51 : '-----------------') 52 53 NFAIL = 0 54 55* PGPLOT initialization. 56 CALL PGBEG (0, '/xwindow', 1, 1) 57 58 NAXISJ = NSPEC 59 CRPIXJ = NAXISJ/2 + 1 60 61 RESTFRQ = 1420.40595D6 62 RESTWAV = C/RESTFRQ 63 X1 = 1D9 64 X2 = 2D9 65 CDELTX = (X2 - X1)/(NAXISJ - 1) 66 CRVALX = X1 + (CRPIXJ - 1.0)*CDELTX 67 WRITE (*, 20) X1*1D-9, X2*1D-9, CDELTX*1D-3 68 20 FORMAT (/,'Linear frequency axis, span:',F4.1,' to',F4.1, 69 : ' (GHz), step:',F8.3,' (kHz)',/,'---------------------', 70 : '-----------------------------------------------------') 71 NFAIL = NFAIL + CLOSURE('WAVE-F2W', 0D0, 0D0, NAXISJ, CRPIXJ, 72 : CDELTX, CRVALX) 73 NFAIL = NFAIL + CLOSURE('VOPT-F2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 74 : CDELTX, CRVALX) 75 NFAIL = NFAIL + CLOSURE('ZOPT-F2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 76 : CDELTX, CRVALX) 77 NFAIL = NFAIL + CLOSURE('AWAV-F2A', 0D0, 0D0, NAXISJ, CRPIXJ, 78 : CDELTX, CRVALX) 79 NFAIL = NFAIL + CLOSURE('VELO-F2V', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 80 : CDELTX, CRVALX) 81 NFAIL = NFAIL + CLOSURE('BETA-F2V', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 82 : CDELTX, CRVALX) 83 84 85 RESTWAV = 700D-9 86 RESTFRQ = C/RESTWAV 87 X1 = 300D-9 88 X2 = 900D-9 89 CDELTX = (X2 - X1)/(NAXISJ - 1) 90 CRVALX = X1 + (CRPIXJ - 1D0)*CDELTX 91 WRITE (*, 30) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9 92 30 FORMAT (/,'Linear vacuum wavelength axis, span:',I4,' to',I4, 93 : ' (nm), step:',F9.6,' (nm)',/,'----------------------', 94 : '----------------------------------------------------') 95 NFAIL = NFAIL + CLOSURE('FREQ-W2F', 0D0, 0D0, NAXISJ, CRPIXJ, 96 : CDELTX, CRVALX) 97 NFAIL = NFAIL + CLOSURE('AFRQ-W2F', 0D0, 0D0, NAXISJ, CRPIXJ, 98 : CDELTX, CRVALX) 99 NFAIL = NFAIL + CLOSURE('ENER-W2F', 0D0, 0D0, NAXISJ, CRPIXJ, 100 : CDELTX, CRVALX) 101 NFAIL = NFAIL + CLOSURE('WAVN-W2F', 0D0, 0D0, NAXISJ, CRPIXJ, 102 : CDELTX, CRVALX) 103 NFAIL = NFAIL + CLOSURE('VRAD-W2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 104 : CDELTX, CRVALX) 105 NFAIL = NFAIL + CLOSURE('AWAV-W2A', 0D0, 0D0, NAXISJ, CRPIXJ, 106 : CDELTX, CRVALX) 107 NFAIL = NFAIL + CLOSURE('VELO-W2V', 0D0, RESTWAV, NAXISJ, CRPIXJ, 108 : CDELTX, CRVALX) 109 NFAIL = NFAIL + CLOSURE('BETA-W2V', 0D0, RESTWAV, NAXISJ, CRPIXJ, 110 : CDELTX, CRVALX) 111 112 113 WRITE (*, 40) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9 114 40 FORMAT (/,'Linear air wavelength axis, span:',I4,' to',I4, 115 : ' (nm), step:',F9.6,' (nm)',/,'----------------------', 116 : '----------------------------------------------------') 117 NFAIL = NFAIL + CLOSURE('FREQ-A2F', 0D0, 0D0, NAXISJ, CRPIXJ, 118 : CDELTX, CRVALX) 119 NFAIL = NFAIL + CLOSURE('AFRQ-A2F', 0D0, 0D0, NAXISJ, CRPIXJ, 120 : CDELTX, CRVALX) 121 NFAIL = NFAIL + CLOSURE('ENER-A2F', 0D0, 0D0, NAXISJ, CRPIXJ, 122 : CDELTX, CRVALX) 123 NFAIL = NFAIL + CLOSURE('WAVN-A2F', 0D0, 0D0, NAXISJ, CRPIXJ, 124 : CDELTX, CRVALX) 125 NFAIL = NFAIL + CLOSURE('VRAD-A2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 126 : CDELTX, CRVALX) 127 NFAIL = NFAIL + CLOSURE('WAVE-A2W', 0D0, 0D0, NAXISJ, CRPIXJ, 128 : CDELTX, CRVALX) 129 NFAIL = NFAIL + CLOSURE('VOPT-A2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 130 : CDELTX, CRVALX) 131 NFAIL = NFAIL + CLOSURE('ZOPT-A2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 132 : CDELTX, CRVALX) 133 NFAIL = NFAIL + CLOSURE('VELO-A2V', 0D0, RESTWAV, NAXISJ, CRPIXJ, 134 : CDELTX, CRVALX) 135 NFAIL = NFAIL + CLOSURE('BETA-A2V', 0D0, RESTWAV, NAXISJ, CRPIXJ, 136 : CDELTX, CRVALX) 137 138 139 RESTFRQ = 1420.40595D6 140 RESTWAV = C/RESTFRQ 141 X1 = -0.96D0*C 142 X2 = 0.96D0*C 143 CDELTX = (X2 - X1)/(NAXISJ - 1) 144 CRVALX = X1 + (CRPIXJ - 1D0)*CDELTX 145 WRITE (*, 50) INT(X1), INT(X2), INT(CDELTX) 146 50 FORMAT (/,'Linear velocity axis, span:',I11,' to',I10, 147 : ' m/s, step:',I6,' (m/s)',/,'-----------------------', 148 : '---------------------------------------------------') 149 NFAIL = NFAIL + CLOSURE('FREQ-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 150 : CDELTX, CRVALX) 151 NFAIL = NFAIL + CLOSURE('AFRQ-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 152 : CDELTX, CRVALX) 153 NFAIL = NFAIL + CLOSURE('ENER-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 154 : CDELTX, CRVALX) 155 NFAIL = NFAIL + CLOSURE('WAVN-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 156 : CDELTX, CRVALX) 157 NFAIL = NFAIL + CLOSURE('VRAD-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 158 : CDELTX, CRVALX) 159 NFAIL = NFAIL + CLOSURE('WAVE-V2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 160 : CDELTX, CRVALX) 161 NFAIL = NFAIL + CLOSURE('VOPT-V2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 162 : CDELTX, CRVALX) 163 NFAIL = NFAIL + CLOSURE('ZOPT-V2W', 0D0, RESTWAV, NAXISJ, CRPIXJ, 164 : CDELTX, CRVALX) 165 NFAIL = NFAIL + CLOSURE('AWAV-V2A', 0D0, RESTWAV, NAXISJ, CRPIXJ, 166 : CDELTX, CRVALX) 167 168 169 RESTWAV = 650D-9 170 RESTFRQ = C/RESTWAV 171 X1 = 300D-9 172 X2 = 1000D-9 173 CDELTX = (X2 - X1)/(NAXISJ - 1) 174 CRVALX = X1 + (CRPIXJ - 1D0)*CDELTX 175 WRITE (*, 60) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9 176 60 FORMAT (/,'Vacuum wavelength grism axis, span:',I4,' to',I5, 177 : ' (nm), step:',F9.6,' (nm)',/,'----------------------', 178 : '----------------------------------------------------') 179 NFAIL = NFAIL + CLOSURE('FREQ-GRI', 0D0, 0D0, NAXISJ, CRPIXJ, 180 : CDELTX, CRVALX) 181 NFAIL = NFAIL + CLOSURE('AFRQ-GRI', 0D0, 0D0, NAXISJ, CRPIXJ, 182 : CDELTX, CRVALX) 183 NFAIL = NFAIL + CLOSURE('ENER-GRI', 0D0, 0D0, NAXISJ, CRPIXJ, 184 : CDELTX, CRVALX) 185 NFAIL = NFAIL + CLOSURE('WAVN-GRI', 0D0, 0D0, NAXISJ, CRPIXJ, 186 : CDELTX, CRVALX) 187 NFAIL = NFAIL + CLOSURE('VRAD-GRI', RESTFRQ, 0D0, NAXISJ, CRPIXJ, 188 : CDELTX, CRVALX) 189 NFAIL = NFAIL + CLOSURE('WAVE-GRI', 0D0, 0D0, NAXISJ, CRPIXJ, 190 : CDELTX, CRVALX) 191 NFAIL = NFAIL + CLOSURE('VOPT-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ, 192 : CDELTX, CRVALX) 193 NFAIL = NFAIL + CLOSURE('ZOPT-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ, 194 : CDELTX, CRVALX) 195 NFAIL = NFAIL + CLOSURE('AWAV-GRI', 0D0, 0D0, NAXISJ, CRPIXJ, 196 : CDELTX, CRVALX) 197 NFAIL = NFAIL + CLOSURE('VELO-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ, 198 : CDELTX, CRVALX) 199 NFAIL = NFAIL + CLOSURE('BETA-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ, 200 : CDELTX, CRVALX) 201 202 203* Reproduce Fig. 5 of Paper III. 204 NAXISJ = 1700 205 CRPIXJ = 719.8D0 206 CRVALX = 7245.2D-10 207 CDELTX = 2.956D-10 208 RESTWAV = 8500D-10 209 RESTFRQ = C/RESTWAV 210 X1 = CRVALX + (1 - CRPIXJ)*CDELTX 211 X2 = CRVALX + (NAXISJ - CRPIXJ)*CDELTX 212 MARS(5) = 0D0 213 MARS(6) = 0D0 214 WRITE (*, 70) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9 215 70 FORMAT (/,'Air wavelength grism axis, span:',I4,' to',I5, 216 : ' (nm), step:',F9.6,' (nm)',/,'----------------------', 217 : '----------------------------------------------------') 218 NFAIL = NFAIL + CLOSURE('AWAV-GRA', 0D0, 0D0, NAXISJ, CRPIXJ, 219 : CDELTX, CRVALX) 220 NFAIL = NFAIL + CLOSURE('VELO-GRA', 0D0, RESTWAV, NAXISJ, CRPIXJ, 221 : CDELTX, CRVALX) 222 223 CALL PGASK(0) 224 CALL PGEND() 225 226 227 IF (NFAIL.NE.0) THEN 228 WRITE (*, 80) NFAIL 229 80 FORMAT (/,'FAIL:',I5,' closure residuals exceed reporting ', 230 : 'tolerance.') 231 ELSE 232 WRITE (*, 90) 233 90 FORMAT (/,'PASS: All closure residuals are within reporting ', 234 : 'tolerance.') 235 END IF 236 237 END 238 239*======================================================================= 240 241 INTEGER FUNCTION CLOSURE (CTYPES, RESTFRQ, RESTWAV, NAXISJ, 242 : CRPIXJ, CDELTX, CRVALX) 243 244 INTEGER NSPEC 245 PARAMETER (NSPEC = 8001) 246 247 INTEGER J, NAXISJ, NFAIL, RESTREQ, STAT1(NSPEC), STAT2(NSPEC), 248 : STATUS 249 REAL TMP, X(NSPEC), XMIN, XMAX, Y(NSPEC), YMAX, YMIN 250 DOUBLE PRECISION CDELTS, CDELTX, CLOS(NSPEC), CRPIXJ, CRVALS, 251 : CRVALX, DSDX, MARS(0:6), RESID, RESIDMAX, RESTFRQ, 252 : RESTWAV, SPEC1(NSPEC), SPEC2(NSPEC), TOL 253 CHARACTER CTYPES*8, PTYPE, SCODE*3, SNAME*21, STYPE*8, TITLE*80, 254 : UNITS*7, XTYPE, YLAB*80 255 256* On some systems, such as Sun Sparc, the struct MUST be aligned 257* on a double precision boundary, done here using an equivalence. 258* Failure to do this may result in mysterious "bus errors". 259 INCLUDE 'spx.inc' 260 INCLUDE 'spc.inc' 261 INTEGER SPC(SPCLEN) 262 DOUBLE PRECISION DUMMY 263 EQUIVALENCE (SPC,DUMMY) 264 265 COMMON /SPECTRO/ MARS 266 267 DATA TOL /1D-11/ 268*----------------------------------------------------------------------- 269 CLOSURE = 0 270 271* Get keyvalues for the required spectral axis type. 272 STATUS = SPCXPS (CTYPES, CRVALX, RESTFRQ, RESTWAV, PTYPE, XTYPE, 273 : RESTREQ, CRVALS, DSDX) 274 IF (STATUS.NE.0) THEN 275 WRITE (*, 10) STATUS, CTYPES 276 10 FORMAT ('ERROR',I2,' from SPCXPS for',A,'.') 277 RETURN 278 END IF 279 CDELTS = CDELTX * DSDX 280 281 STATUS = SPCINI(SPC) 282 283 IF (CTYPES(6:6).EQ.'G') THEN 284* KPNO MARS spectrograph grism parameters. 285 STATUS = SPCPTD (SPC, SPC_PV, MARS(0), 0) 286 STATUS = SPCPTD (SPC, SPC_PV, MARS(1), 1) 287 STATUS = SPCPTD (SPC, SPC_PV, MARS(2), 2) 288 STATUS = SPCPTD (SPC, SPC_PV, MARS(3), 3) 289 STATUS = SPCPTD (SPC, SPC_PV, MARS(4), 4) 290 STATUS = SPCPTD (SPC, SPC_PV, MARS(5), 5) 291 STATUS = SPCPTD (SPC, SPC_PV, MARS(6), 6) 292 END IF 293 294* Construct the axis. 295 DO 20 J = 1, NAXISJ 296 SPEC1(J) = (J - CRPIXJ)*CDELTS 297 20 CONTINUE 298 299 WRITE (*, 30) CTYPES, CRVALS+SPEC1(1), CRVALS+SPEC1(NAXISJ), 300 : CDELTS 301 30 FORMAT (A,' (CRVALk+w) range: ',1PE13.6,' to ',1PE13.6,', step: ', 302 : 1PE13.6) 303 304 305* Initialize. 306 STATUS = SPCPTI (SPC, SPC_FLAG, 0, 0) 307 STATUS = SPCPTD (SPC, SPC_CRVAL, CRVALS, 0) 308 STATUS = SPCPTD (SPC, SPC_RESTFRQ, RESTFRQ, 0) 309 STATUS = SPCPTD (SPC, SPC_RESTWAV, RESTWAV, 0) 310 STATUS = SPCPTC (SPC, SPC_TYPE, CTYPES, 0) 311 STATUS = SPCPTC (SPC, SPC_CODE, CTYPES(6:8), 0) 312 313* Convert the first to the second. 314 STATUS = SPCX2S(SPC, NAXISJ, 1, 1, SPEC1, SPEC2, STAT1) 315 IF (STATUS.NE.0) THEN 316 WRITE (*, 40) STATUS 317 40 FORMAT ('SPCX2S ERROR',I2,'.') 318 RETURN 319 END IF 320 321* Convert the second back to the first. 322 STATUS = SPCS2X(SPC, NAXISJ, 1, 1, SPEC2, CLOS, STAT2) 323 IF (STATUS.NE.0) THEN 324 WRITE (*, 50) STATUS 325 50 FORMAT ('SPCS2X ERROR',I2,'.') 326 RETURN 327 END IF 328 329 330* Test closure. 331 NFAIL = 0 332 RESIDMAX = 0D0 333 STATUS = SPCGTC (SPC, SPC_TYPE, STYPE, 0) 334 DO 90 J = 1, NAXISJ 335 IF (STAT1(J).NE.0) THEN 336 WRITE (*, 60) CTYPES, SPEC1(J), STYPE, STAT1(J) 337 60 FORMAT (A,': w =',1PE20.12,' -> ',A4,' = ???, stat = ',I2) 338 GO TO 90 339 END IF 340 341 IF (STAT2(J).NE.0) THEN 342 WRITE (*, 70) CTYPES, SPEC1(J), STYPE, SPEC2(J), STAT2(J) 343 70 FORMAT (A,': w =',1PE20.12,' -> ',A4,' =',1PE20.12, 344 : ' -> w = ???, stat = ',I2) 345 GO TO 90 346 END IF 347 348 RESID = ABS((CLOS(J) - SPEC1(J))/CDELTS) 349 IF (RESID.GT.RESIDMAX) RESIDMAX = RESID 350 351 IF (RESID.GT.TOL) THEN 352 NFAIL = NFAIL + 1 353 WRITE (*, 80) CTYPES, SPEC1(J), STYPE, SPEC2(J), CLOS(J), 354 : RESID 355 80 FORMAT (A,': w =',1PE20.12,' -> ',A4,' =',1PE20.12,' ->',/, 356 : ' w =',1PE20.12,', resid =',1PE8.1) 357 END IF 358 90 CONTINUE 359 360 WRITE (*, 100) CTYPES, RESIDMAX 361 100 FORMAT (A,': Maximum closure residual =',1PE8.1,' pixel') 362 363 364* Draw graph. 365 CALL PGBBUF() 366 CALL PGERAS() 367 368 XMIN = REAL(CRVALS + SPEC1(1)) 369 XMAX = REAL(CRVALS + SPEC1(NAXISJ)) 370 YMIN = REAL(SPEC2(1)) - XMIN 371 YMAX = YMIN 372 DO 110 J = 1, NAXISJ 373 X(J) = REAL(J) 374 Y(J) = REAL(SPEC2(J) - (CRVALS + SPEC1(J))) 375 IF (Y(J).GT.YMAX) YMAX = Y(J) 376 IF (Y(J).LT.YMIN) YMIN = Y(J) 377 110 CONTINUE 378 379 J = INT(CRPIXJ+1) 380 IF (Y(J).LT.0D0) then 381 TMP = YMIN 382 YMIN = YMAX 383 YMAX = TMP 384 END IF 385 386 CALL PGASK(0) 387 CALL PGENV(1.0, REAL(NAXISJ), YMIN, YMAX, 0, -1) 388 389 CALL PGSCI(1) 390 CALL PGBOX('ABNTS', 0.0, 0, 'BNTS', 0.0, 0) 391 392 STATUS = SPCTYP (CTYPES, STYPE, SCODE, SNAME, UNITS, PTYPE, XTYPE, 393 : RESTREQ) 394 DO 120 J = 21, 1, -1 395 IF (SNAME(J:J).NE.' ') GO TO 130 396 120 CONTINUE 397 130 YLAB = SNAME(:J) // ' - correction ' // UNITS 398 TITLE = CTYPES // ': CRVALk + w ' // UNITS 399 CALL PGLAB('Pixel coordinate', YLAB, TITLE) 400 401 CALL PGAXIS('N', 0.0, YMAX, REAL(NAXISJ), YMAX, XMIN, XMAX, 0.0, 402 : 0, -0.5, 0.0, 0.5, -0.5, 0.0) 403 404 CALL PGAXIS('N', REAL(NAXISJ), ymin, REAL(NAXISJ), YMAX, 405 : REAL(YMIN/CDELTS), REAL(YMAX/CDELTS), 0.0, 0, 0.5, 406 : 0.0, 0.5, 0.1, 0.0) 407 CALL PGMTXT('R', 2.2, 0.5, 0.5, 'Pixel offset') 408 409 CALL PGLINE(NAXISJ, X, Y) 410 CALL PGSCI(7) 411 CALL PGPT1(REAL(CRPIXJ), 0.0, 24) 412 CALL PGEBUF() 413 414 WRITE (*, '(A,$)') 'Type <RETURN> for next page: ' 415 READ (*, *, END=140) 416 140 WRITE (*, *) 417 418 CLOSURE = NFAIL 419 420 END 421