forked from susanasanchez/GUIPSY
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrfits.c
2841 lines (2615 loc) · 113 KB
/
rfits.c
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
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* rfits.c
Copyright (c) Kapteyn Laboratorium Groningen 1990
All Rights Reserved.
#> rfits.dc1
Program: RFITS
Purpose: Load FITS files from tape into a GIPSY set.
Category: FITS, TAPES, UTILITY
File: rfits.c
Author: K.G. Begeman
Description: The program RFITS is designed to load FITS files into
GIPSY sets. Since the F in FITS stands for flexible,
it is not that simple to comprehend the data structure
which is to be loaded. Therefore RFITS needs a lot of
interaction with the user, and the user must know what
is on tape and where it is on tape.
Since this program must also be flexible, it will always
be in devellopment, which means that there will be
a constant interaction between the user and the programmer.
So don't hesitate to bring forward suggestions for
improvement, or even criticism.
Keywords:
*** AUTO= Try to run in automatice mode [Y].
In automatice mode the program will do its best to
comprehend the FITS structure on tape. FITS items for
which the program has some reasonable defaults will
not be asked. In non-automatic mode all items will be
asked.
*** NODATA= Create only set headers [N].
With NODATA=Y no data is read, only the set headers
are created.
*** HISTORY= Copy history information from FITS file to Set header [Y].
*** COMMENT= Copy comments from FITS file to Set Header [Y].
FITSFILE= Name of file with FITS data [read from tape device].
If a file name is entered, the keywords INTAPE= and
INFILES= will be skipped. This keyword, and all the
following keywords are repeated. After the first filename
is read, the default will change to [RFITS quits].
INTAPE= Tape device to load from [list of all tape devices].
INFILES= File numbers on tape to load data from [RFITS quits].
Maximum number of files is 200. The first time this
keyword is asked, the current file on tape has file
number 0, the previous file (if present) has file number
-1 and the next file has file number +1. All numbers will
be relative to this initial position. It is also important
to know that the program reads the tape in forward direction,
meaning that it does not care about the order of file numbers
entered by the user.
This keyword, and all the following keywords are repeated.
OUTSET= Set to load data into. The following keywords depend
on whether the set already exists or not.
*** BLANK= Redefine the FITS BLANK value as specified in the FITS
header [value from FITS header]. This option is useful
if the program which created the FITS file made a
mistake.
--------------------------------------------------------
The following keywords deal with the FITS structure on
tape. If the program comprehends this structure, the
keywords will be hidden. You can force all keywords to
be prompted with AUTO=N.
The asterix in the keyword will be replaced by the
corresponding axis number.
--------------------------------------------------------
CTYPE*= Enter the corresponding axis name. Default is chosen as
smart as possible. If the set does not exist, the user may
type any other legal axis name. If the axis name found in
the FITS file is NOT a correct axis name, the user must
supply a legal axis name.
It is also possible to SKIP an axis or to HIDE an axis.
The default will again be chosen as smart as possible.
NOTE: Frequency axes to which a velocity system is
attached must be indicated by FREQ-OHEL, FREQ-OLSR,
FREQ-RHEL or FREQ-RLSR.
NOTE: This keyword is also asked (hidden) if the set
does not yet exist. It will ask for extra axes names
to be added to the new set, all with length 1.
CUNIT*= Enter the units of the axis on the FITS tape. This
keyword only appears if the set does not yet exist
and the corresponding FITS descriptor is not present.
Default depends on CTYPE and is normally as prescribed
by the FITS standard.
CRVAL*= Enter the reference value on the axis. This keyword only
appears if the set does not yet exist.
Default will be chosen as best we can.
CDELT*= Enter the pixel separation along the axis. This keyword only
appears if the set does not yet exist.
Note: WSRT data might have the wrong sign for CDELT3.
*** LIMITS*= Enter the range in grids along the axis [whole axis].
If the set does already exist, this keyword will not be asked.
Note that GIPSY does not allow to extend axes in any
direction. Only the last axis in a GIPSY set may be
extended, so with this keyword you can make room for
data which should be loaded at lower grid positions.
GRID*= Enter the first grid position along the axis [first on axis].
This keyword is only asked if the set does not yet exist.
Default will be chosen as best we can.
FREQ0= Rest frequency of observation in Hz. This keyword only
appears if FITS descriptor not found in FITS file and
the axis is a frequency axis with velocity information.
VEL= Velocity of first grid on axis. This keywords is only asked
for old FITS files which have been written by VMS GIPSY.
Example: <USER > RFITS,AUTO=N
<USER > RFITS,INTAPE=IMPORT
<USER > RFITS,INFILES=31
<USER > RFITS,OUTSET=NGC4214
<USER > RFITS,CTYPE1=RA-NCP
<USER > RFITS,CUNIT1=
<USER > RFITS,LIMITS1=
<USER > RFITS,GRID1=
<USER > RFITS,CTYPE2=DEC-NCP
<USER > RFITS,CUNIT2=
<USER > RFITS,LIMITS2=
<USER > RFITS,GRID2=
<USER > RFITS,CTYPE3=FREQ-OHEL
<USER > RFITS,CUNIT3=MHZ
<USER > RFITS,CDELT3=-2.5/64
<USER > RFITS,LIMITS3=1 32
<USER > RFITS,GRID3=32
<USER > RFITS,FREQ0=
<USER > RFITS,CTYPE4=
<USER > RFITS,CTYPE5=
<USER > RFITS,INFILES=
Updates: May 24, 1990: KGB, Document created.
Nov 26, 1991: KGB, Automatic mode implemented.
Feb 26, 1993: KGB, Modification for NEWSTAR FITS.
Mar 19, 1993: KGB, Extending the dimensions of new set.
Apr 23, 1993: KGB, Bug repaired, keyword NODATA added.
Apr 14, 1994: KGB, Bug when reading multiple NMAP cubes repaired.
Feb 5, 1996: KGB, Keyword FITSFILE= implemented.
Dec 14, 1999: VOG, Allow axis names X1...Xn as default for
missing CTYPEs in FITS header
Mar 3, 2005: JPT, Implemented 64-bit float (BITPIX = -64)
Apr 12, 2009: VOG, Replaced definition of function nint with
one that uses function floor(). This is to
enforce compatibility with other coordinate
routines and to avoid problems with values
of CRPIX that end on .5
#<
*/
/*
* include files:
*/
#include "ctype.h" /* <ctype.h> */
#include "float.h" /* <float.h> */
#include "limits.h" /* <limits.h> */
#include "math.h" /* <math.h> */
#include "stdio.h" /* <stdio.h> */
#include "stdlib.h" /* <stdlib.h> */
#include "string.h" /* <string.h> */
#include "gipsyc.h" /* GIPSY definitions and symbols */
#include "cmain.h" /* C MAIN PROGRAM */
#include "anyout.h" /* define anyout_c */
#include "axtype.h" /* define axtype_c */
#include "cancel.h" /* define cancel_c */
#include "clspfp.h" /* define clspfp_c */
#include "dcdint.h" /* define dcdint_c */
#include "error.h" /* define error_c */
#include "factor.h" /* define factor_c */
#include "finis.h" /* define finis_c */
#include "ftsd_find.h" /* define ftsd_find_c */
#include "ftsd_rchar.h" /* define ftsd_rchar_c */
#include "ftsd_rdble.h" /* define ftsd_rdble_c */
#include "ftsd_rint.h" /* define ftsd_rint_c */
#include "ftsd_rlog.h" /* define ftsd_rlog_c */
#include "gds_close.h" /* define gds_close_c */
#include "gds_create.h" /* define gds_create_c */
#include "gds_exist.h" /* define gds_exist_c */
#include "gds_extend.h" /* define gds_extend_c */
#include "gdsc_grid.h" /* define gdsc_grid_c */
#include "gdsc_substruct.h" /* define gdsc_substruct_c */
#include "gdsc_word.h" /* define gdsc_word_c */
#include "gdsd_length.h" /* define gdsd_length_c */
#include "gdsd_rchar.h" /* define gdsd_rchar_c */
#include "gdsd_rdble.h" /* define gdsd_rdble_c */
#include "gdsd_rint.h" /* define gdsd_rint_c */
#include "gdsd_wchar.h" /* define gdsd_wchar_c */
#include "gdsd_wdble.h" /* define gdsd_wdble_c */
#include "gdsd_wfits.h" /* define gdsd_wfits_c */
#include "gdsi_write.h" /* define gdsi_write_c */
#include "gdsd_wvar.h" /* define gdsd_wvar_c */
#include "init.h" /* define init_c */
#include "initptr.h" /* define initptr_c */
#include "insideptr.h" /* define insideptr_c */
#include "mtbsf.h" /* define mtbsf_c */
#include "mtbsr.h" /* define mtbsr_c */
#include "mtclose.h" /* define mtclose_c */
#include "mtfsf.h" /* define mtfsf_c */
#include "mtname.h" /* define mtname_c */
#include "mtopen.h" /* define mtopen_c */
#include "mtreadc.h" /* define mtreadc_c */
#include "mtrew.h" /* define mtrew_c */
#include "nelc.h" /* define nelc_c */
#include "reject.h" /* define reject_c */
#include "setfblank.h" /* define setfblank_c */
#include "sortia.h" /* define sortia_c */
#include "spfpfl.h" /* define spfpfl_c */
#include "dpfpfl.h" /* define dpfpfl_c */
#include "status.h" /* define status_c */
#include "usercharu.h" /* define usercharu_c */
#include "userdble.h" /* define userdble_c */
#include "userint.h" /* define userint_c */
#include "userlog.h" /* define userlog_c */
#include "usertext.h" /* define usertext_c */
#include "velpro.h" /* define velpro_c */
/*
* Define the keywords and associated messages.
*/
#define VERSION "1.3" /* change version number on this line */
#define AUTO_KEY tofchar("AUTO=")
#define AUTO_MES tofchar("Automatic mode [Y]")
#define COMMENT_KEY tofchar("COMMENT=")
#define COMMENT_MES tofchar("Copy FITS comments? [Y]")
#define FITSFILE_KEY tofchar("FITSFILE=")
#define FITSFILE_MES1 tofchar("Name of FITS file [read from tape device]")
#define FITSFILE_MES2 tofchar("Name of FITS file [program quits]")
#define FMAKE_ERR tofchar("Cannot allocate enough memory!")
#define HISTORY_KEY tofchar("HISTORY=")
#define HISTORY_MES tofchar("Copy FITS history? [Y]")
#define INTAPE_DEV tofchar("?INTAPE=Tape device to load from [list of all tape devices]")
#define INTAPE_ERR1 tofchar("No such tape device!")
#define INFILES_KEY tofchar("INFILES=")
#define INFILES_MES tofchar("Give file numbers on tape [program quits]")
#define MAIN_ERR1 tofchar("Program expects 8 bit bytes!")
#define MAIN_ERR2 tofchar("Program expects 1 byte chars!")
#define MAIN_ERR3 tofchar("Program expects 2 byte shorts!")
#define MAIN_ERR4 tofchar("Program expects 4 byte ints!")
#define MAIN_ERR5 tofchar("Program expects 4 byte floats!")
#define NODATA_KEY tofchar("NODATA=")
#define NODATA_MES tofchar("Skip reading in the image part? [N]")
#define OKAY_KEY tofchar("OKAY=")
#define OKAY_MES tofchar("Set already exists. Overwrite? [Y]/N")
#define OUTSET_KEY tofchar("OUTSET=")
#define OUTSET_MES tofchar("Set to load data into")
#define POSITION_ERR tofchar("error positioning tape!")
#define READHEADER_ERR1 tofchar("End Of Information encountered!")
#define READHEADER_ERR2 tofchar("Wrong block size!")
#define READHEADER_ERR3 tofchar("Tape error ?!?")
#define READHEADER_ERR4 tofchar("Error reading FITS header!")
#define DECODEHEAD_ERR1 tofchar("FITS descriptor SIMPLE not present!")
#define DECODEHEAD_ERR2 tofchar("SIMPLE = FALSE (will try anyway!)")
#define DECODEHEAD_ERR3 tofchar("FITS descriptor BITPIX nor present!")
#define DECODEHEAD_ERR4 tofchar("BITPIX not compatible with this machine!")
#define DECODEHEAD_ERR5 tofchar("Descriptor NAXIS not present!")
#define GAMBLE1_ERR1 tofchar("Cannot hide this axis! Try again!")
#define GAMBLE1_ERR2 tofchar("Cannot skip this axis! Try again!")
#define GAMBLE1_ERR3 tofchar("Unknown axis! Try again!")
#define GAMBLE1_ERR4 tofchar("Illegal limits! Try again!")
#define GAMBLE1_ERR5 tofchar("No data in range! Try again!")
#define GAMBLE2_ERR1 tofchar("Axis not in set!")
#define GAMBLE2_ERR2 tofchar("GRID outside range!")
#define GAMBLE2_ERR3 tofchar("Error extending set!")
#define GAMBLE2_ERR4 tofchar("Illegal limits! Try again!")
#define GAMBLE2_ERR5 tofchar("No data in range! Try again!")
#define COPYDATA_ERR1 tofchar("GRID outside range!" )
#define COPYDATA_ERR2 tofchar("Error extending set!")
#define COPYDATA_ERR3 tofchar("BITPIX not compatible with this machine!")
#define COPYDATA_ERR4 tofchar("End Of Information encountered!")
#define COPYDATA_ERR5 tofchar("Wrong block size!")
#define COPYDATA_ERR6 tofchar("Tape error ?!?")
#define COPYDATA_ERR7 tofchar("Error skipping to next file!")
/*
* define some parameters and functions:
*/
#define MUNIT 1 /* mask for CUNIT/DUNIT */
#define MTYPE 2 /* mask for CTYPE/DTYPE */
#define MRVAL 4 /* mask for CRVAL/DRVAL */
#define MRPIX 8 /* mask for CRPIX/DRPIX */
#define MROTA 16 /* mask for CROTA/DROTA */
#define MDELT 32 /* mask for CDELT/DDELT */
#define DSC_SKIP 1 /* skip this descriptor */
#define DSC_SUB 2 /* descriptor at subset level */
#define DSC_CHECK 4 /* dont't overwrite old descriptor */
#define DSC_COMMENT 8 /* comment */
#define DSC_HISTORY 16 /* history */
#define MAXAXES 16 /* maximum number of axes possible */
/* size of descriptor list */
#define DSCLISTLEN (sizeof(dsc_buf)/sizeof(dsc_struct))
#define MAXDEVNAMLEN (FILENAME_MAX) /* maximum length of tape name */
#define MAXFTSDSCLEN 8 /* length of FITS descriptor */
#define MAXFTSNAMLEN 18 /* length of FITS character value */
#define MAXFTSCHRLEN 78 /* .. */
#define MAXSETNAMLEN 80 /* maximum length of set name */
#define MAXSYNONYMS 100 /* maximum numver of synonyms */
//#define MAXTEXTLEN 128 /* maximum length of text strings */
#define MAXTEXTLEN 512 /* maximum length of text strings */
#define MAXFILES 200 /* maximum number of files */
#define FITSRECLEN 80 /* length of a FITS record on tape */
#define FITSBLOCKLEN 2880 /* length of FITS blocks on tape */
#define FITSBLOCKFACTOR 1 /* FITS blocking factor (STANDARD) */
/*
* some macros:
*/
#define READ( r, b ) \
if ( fid != NULL ) { \
tst = fread( b.a, sizeof( char ), b.l, fid ); \
} else { \
tst = mtreadc_c( &mtid, b ); \
}
#define PANIC( s ) \
{ \
fint elev = 4; \
if ( fid == NULL ) { \
(void) mtrew_c( &mtid ); \
(void) mtclose_c( &mtid ); \
} else { \
(void) fclose( fid ); \
} \
error_c( &elev, s ); \
}
/*
* define the structures:
*/
typedef struct { /* structure for axis info */
fint naxis; /* length of an axis */
fint pmask; /* mask for primary axis */
double cdelt; /* increment in units along axis */
double crota; /* rotation angle of axis */
double crpix; /* reference pixel of axis */
double crval; /* reference value at reference pixel */
char ctypeb[MAXFTSNAMLEN]; /* axis name */
char cunitb[MAXFTSNAMLEN]; /* units of axis */
fchar ctype; /* f character points to ctypeb */
fchar cunit; /* f character points to cunitb */
fint smask; /* mask for secondary axis */
double ddelt; /* increment along secondary axis */
double drota; /* secondary rotation angle of axis */
double drpix; /* secondary reference pixel of axis */
double drval; /* units at secondary reference pixel */
char dtypeb[MAXFTSNAMLEN]; /* secondary axis name */
char dunitb[MAXFTSNAMLEN]; /* secondary units of axis */
fchar dtype; /* f character points to dtypeb */
fchar dunit; /* f character points to dunitb */
fint matchax; /* matching axis */
fint axtype; /* type of axis */
fint skysys; /* sky system code */
fint prosys; /* projection system code */
fint velsys; /* velocity system code */
fint low; /* lower grid value on axis */
fint upp; /* upper grid value on axis */
fint min; /* wanted lower grid value on axis */
fint max; /* wanted upper grid value on axis */
fint grid; /* grid position */
} ax_struct;
typedef struct { /* set info structure */
ax_struct ax[MAXAXES]; /* contains axis info */
char setb[MAXSETNAMLEN]; /* name of set */
char instrumeb[MAXFTSNAMLEN]; /* name of instrument */
fchar set; /* f character points to setb */
fchar instrume; /* f character points to instrumeb */
fint subset; /* level of subset */
double epoch; /* epoch of set */
double freq0; /* rest frequency of set */
fint exist; /* does set exist */
fint naxis; /* number of real axes */
fint maxis; /* total number of axes */
float blank; /* blank value */
} set_struct;
typedef struct { /* fits info structure */
ax_struct ax[MAXAXES]; /* contains axis info */
char devb[MAXDEVNAMLEN+1]; /* buffer for tape device name */
fchar dev; /* f character points to devb */
char instrumeb[MAXFTSNAMLEN]; /* name of instrument */
fchar instrume; /* f character points to instrumeb */
double bscal; /* data scaling factor */
double bzero; /* data offset factor */
double epoch; /* epoch of set */
double freq0; /* rest frequency of set */
fint axused; /* number of axes in used already */
fint nbytes; /* number of bytes / pixel */
fint bitpix; /* number of bits / pixel */
fint naxis; /* number of real axes */
fint maxis; /* total number of axes */
fint blank; /* blank value on tape */
} fts_struct;
typedef struct { /* synonym structure */
char cnameb[MAXFTSNAMLEN]; /* buffer for synonym */
fchar cname; /* points to cnameb */
char ctypeb[MAXFTSNAMLEN]; /* buffer for ctype */
fchar ctype; /* points to ctypeb */
} syn_struct; /* THE STRUCT */
typedef struct { /* descriptor structure */
char *dsc; /* descriptor name */
char *rep; /* replace by ... */
int code; /* code */
} dsc_struct; /* THE STRUCT */
typedef unsigned char byte; /* make'em unsigned */
/*
* define the global variables:
*/
static fts_struct fts_buf; /* contains FITS info */
static set_struct set_buf; /* contains GDS info */
static syn_struct syn_buf[MAXSYNONYMS]; /* contains synonyms */
static dsc_struct dsc_buf[] = {
{ "APLAB" , NULL , DSC_SUB }, /* tape label of Antenna Pattern */
{ "APSET" , NULL , DSC_SKIP }, /* set of AP */
{ "APSSET" , NULL , DSC_SKIP }, /* subset of AP */
{ "APVSN" , NULL , DSC_SUB }, /* tape volume of Antenna Pattern */
{ "BANDW" , NULL , DSC_CHECK }, /* total bandwidth of observation */
{ "BITPIX" , NULL , DSC_SKIP }, /* for encoding only */
{ "BLANK" , NULL , DSC_SKIP }, /* blank value */
{ "BLGRAD" , NULL , DSC_CHECK }, /* baseline grading function */
{ "BLOCKED" , NULL , DSC_SKIP }, /* only for tapes */
{ "BMMAJ" , NULL , DSC_CHECK }, /* major axis of beam */
{ "BMMIN" , NULL , DSC_CHECK }, /* minor axis of beam */
{ "BMPA" , NULL , DSC_CHECK }, /* PA major axis beam */
{ "BSCALE" , NULL , DSC_SKIP }, /* scaling factor */
{ "BUNIT" , NULL , DSC_CHECK }, /* data units (WU,MJY/SR,...) */
{ "BZERO" , NULL , DSC_SKIP }, /* zero level */
{ "CDELT" , NULL , DSC_SKIP }, /* primary grid separation */
{ "COMMENT" , NULL , DSC_COMMENT | DSC_SUB }, /* Comment records */
{ "CROTA" , NULL , DSC_SKIP }, /* primary rotation angle */
{ "CRPIX" , NULL , DSC_SKIP }, /* primary reference pixel */
{ "CRVAL" , NULL , DSC_SKIP }, /* primary reference value */
{ "CTYPE" , NULL , DSC_SKIP }, /* primary axis name */
{ "CUNIT" , NULL , DSC_SKIP }, /* primary axis units */
{ "DATAMAX" , NULL , DSC_SUB }, /* maximum in FITS file */
{ "DATAMIN" , NULL , DSC_SUB }, /* minimum in FITS file */
{ "DATE" , NULL , DSC_CHECK }, /* tape writing date */
{ "DATE-OBS", NULL , DSC_CHECK }, /* date of observation */
{ "DDELT" , NULL , DSC_SKIP }, /* secondary grid separation */
{ "DROTA" , NULL , DSC_SKIP }, /* secondary rotation angle */
{ "DRPIX" , NULL , DSC_SKIP }, /* secondary reference pixel */
{ "DRVAL" , NULL , DSC_SKIP }, /* secondary referenc value */
{ "DTYPE" , NULL , DSC_SKIP }, /* secondary axis name */
{ "DUNIT" , NULL , DSC_SKIP }, /* secondary axis units */
{ "EPOCH" , NULL , DSC_SKIP }, /* epoche of observation (years) */
{ "EQUINOX" , NULL , DSC_SKIP }, /* equinox */
{ "EXTEND" , NULL , DSC_SKIP }, /* only for tapes */
{ "FILEID" , NULL , DSC_SKIP }, /* identification of file */
{ "FREQ0" , NULL , DSC_SKIP }, /* Rest frequency of spectral line (Hz) */
{ "FSCDEC" , NULL , DSC_CHECK }, /* DEC fringe stopping center (degrees) */
{ "FSCRA" , NULL , DSC_CHECK }, /* RA fringe stopping center (degrees) */
{ "GRIDTYPE", NULL , DSC_CHECK }, /* type of grid */
{ "HISTORY" , NULL , DSC_HISTORY | DSC_SUB }, /* History records */
{ "INSTRUME", NULL , DSC_SKIP }, /* name of instrument */
{ "LMAX" , NULL , DSC_SKIP }, /* ????? */
{ "LMIN" , NULL , DSC_SKIP }, /* ????? */
{ "MAPCODE" , NULL , DSC_SKIP }, /* LINEMAP map code */
{ "MAPLAB" , NULL , DSC_SUB }, /* tape label of map archive */
{ "MAPVSN" , NULL , DSC_SUB }, /* tape volume of map archive */
{ "MAXBASE" , NULL , DSC_CHECK }, /* maximum baseline */
{ "MAXINT" , NULL , DSC_SKIP }, /* ????? */
{ "MINBASE" , NULL , DSC_CHECK }, /* minimum baseline */
{ "MININT" , NULL , DSC_SKIP }, /* ????? */
{ "MMAX" , NULL , DSC_SKIP }, /* ????? */
{ "MMIN" , NULL , DSC_SKIP }, /* ????? */
{ "NAXIS" , NULL , DSC_SKIP }, /* number of axes */
{ "NBLANK" , NULL , DSC_SUB }, /* number of undefined values in FITS file */
{ "NINTF" , NULL , DSC_CHECK }, /* number of interferometers used */
{ "NOISE" , NULL , DSC_SUB }, /* r.m.s. noise in FITS file */
{ "NORM" , NULL , DSC_SUB }, /* normalizing factor in FFT */
{ "NPOL" , NULL , DSC_CHECK }, /* number of polarizations used */
{ "NRFREQ" , NULL , DSC_CHECK }, /* number of frequency points used */
{ "OBJECT" , NULL , DSC_CHECK }, /* name of object */
{ "OBSDEC" , "PCDEC" , DSC_CHECK }, /* replace with PCDEC */
{ "OBSERVER", NULL , DSC_CHECK }, /* name of observer */
{ "OBSRA" , "PCRA" , DSC_CHECK }, /* replace with PCRA */
{ "OBSTIME" , NULL , DSC_CHECK }, /* observation time */
{ "OBSTYPE" , NULL , DSC_CHECK }, /* type of observation (LINE,CONT) */
{ "ORIGIN" , NULL , DSC_CHECK }, /* tape writing institute */
{ "PCDEC" , NULL , DSC_CHECK }, /* pointing centre DEC */
{ "PCRA" , NULL , DSC_CHECK }, /* pointing centre RA */
{ "REDCODE" , NULL , DSC_SKIP }, /* LINEMAP reduction code */
{ "RESOL" , NULL , DSC_SKIP }, /* spectral resolution */
{ "RESTFREQ", NULL , DSC_SKIP }, /* is in freq0 */
{ "SIMPLE" , NULL , DSC_SKIP }, /* simple must be true */
{ "TAPER" , NULL , DSC_CHECK }, /* frequency taper */
{ "TELESCOP", NULL , DSC_SKIP }, /* name of telescope */
{ "UVBANDW" , NULL , DSC_SKIP }, /* UV bandwidth */
{ "UVFREQ" , NULL , DSC_SKIP }, /* UV frequency */
{ "UVGRID" , NULL , DSC_CHECK }, /* convolving function code */
{ "VELCODE" , NULL , DSC_SKIP }, /* velocity system */
{ "VELREF" , NULL , DSC_SKIP }, /* aips velref */
};
static FILE *fid = NULL; /* file id */
static bool amode = TRUE; /* automatic mode */
static bool comment = TRUE; /* save comment */
static bool history = TRUE; /* save history */
static bool nodata = FALSE; /* load image data */
static fchar fts_head = { NULL, 0 }; /* dynamic memory for FITS header */
static fchar fts_data = { NULL, 0 }; /* dynamic memory for FITS image */
static fchar fts_tape = { NULL, 0 }; /* dynamic memory for FITS tape io */
static fint currentfile = 0; /* current file number */
static fint files[MAXFILES]; /* file numbers */
static fint indata; /* number of bytes in data buffer */
static fint inhead; /* number of bytes in header buffer */
static fint mtid; /* id of tape device */
static fint nfile; /* number of files to load from tape */
static fint nsyns = 0; /* number of synonyms in syn_buf */
static fint nint( double x )
{
/* Pre Apr 2009 def.: return(x > 0.0 ? (fint) ( x + 0.5 ) : (fint) ( x - 0.5 )); */
return( (fint) floor(x + 0.5) );
}
static void fcopy( fchar dest, fchar source)
/*
* fcopy copies a f character from source to dest, filling dest with
* trailing blanks if dest longer than source, or truncating dest if
* dest shorter than source.
*/
{
int n; /* counter */
/* copy loop */
for (n = 0; n < dest.l && n < source.l; n++) {
dest.a[n] = source.a[n]; /* copy character */
}
while (n < dest.l) dest.a[n++] = ' '; /* trailing blanks */
}
static int fcomp( fchar a, fchar b)
/*
* fcomp compares f characters similar to strcmp.
*/
{
return( strncmp( a.a, b.a, a.l > b.l ? b.l : a.l ) );
}
static fchar finit( char *b, int len )
/*
* finit returns a f character which points to the character buffer
* b, which contains len bytes.
* The buffer b is filled with blanks.
*/
{
fchar r; /* return value */
int i; /* counter */
r.a = b; /* assign buffer */
r.l = len; /* length of buffer */
for (i = 0; i < r.l; r.a[i++] = ' '); /* fill with blanks */
return( r ); /* retrun f character */
}
static fchar fmake( fchar f, fint len )
/*
* fmake dynamically allocates len bytes of memory for the fortran
* character f. This is only done if len greater than the current size
* of f.
*/
{
fchar r; /* return value */
if (f.l < len) { /* enlarge */
r.a = realloc( f.a, len ); /* allocate len bytes */
if (r.a == NULL) {
fint elev = 4; /* error code (FATAL) */
error_c( &elev, FMAKE_ERR ); /* error message */
}
r.l = len; /* new length of f */
} else {
r.a = f.a; r.l = f.l; /* use old f character */
}
return( r ); /* return to caller */
}
static void fixup( fchar fc )
/*
* fixup scans the string fc for non_printable characters and leading
* blanks.
*/
{
int k;
for ( k = 0; k < fc.l; k++ ) {
if (!isprint( fc.a[k] )) fc.a[k] = ' ';
}
for ( k = 0; k < fc.l && fc.a[k] == ' '; k++ );
if (k) {
int m = 0;
while ( k < fc.l ) fc.a[m++] = fc.a[k++];
while ( m < fc.l ) fc.a[m++] = ' ';
}
}
static fchar descr( char *name, fint axnum )
/*
* descr creates a fits descriptor consisting of name plus the axis
* number attached. The address of the f character points to a static
* area which is overwritten each call. descr does NOT check whether
* the resulting descriptor name will overflow the static buffer area!
*/
{
static char dscbuf[MAXFTSDSCLEN+1]; /* buffer for descriptor */
fchar r; /* return value */
fint l; /* counter */
l = sprintf( dscbuf, "%s%d", name, axnum ); /* encode descriptor */
while (l < MAXFTSDSCLEN) dscbuf[l++] = ' '; /* blank fill */
r.a = dscbuf; r.l = l; /* initialize return value */
return( r ); /* return to caller */
}
static fint matchaxis( fchar ctype )
/*
* matchaxis compares ctype with the axis names in the existing set.
* For hidden axis, matchaxis returns -2, for other matching axes
* the axtype code is returned, If no matching axis, zero is returned.
*/
{
char cunitb[MAXFTSNAMLEN]; /* buffer for cunit */
char dunitb[MAXFTSNAMLEN]; /* buffer for dunit */
fchar cunit; /* points to cunitb */
fchar dunit; /* points to dunitb */
fint axmatch = 0; /* reset matching axis number */
fint axtype; /* axis type of FITS axis */
fint n; /* loop counter */
fint prosys; /* projection system */
fint skysys; /* sky system */
fint velsys; /* velocity system */
cunit = finit( cunitb, MAXFTSNAMLEN ); /* initialize cunit */
dunit = finit( dunitb, MAXFTSNAMLEN ); /* initialize dunit */
/* obtain axis type */
axtype = axtype_c( ctype, cunit, dunit, &skysys, &prosys, &velsys );
/* search loop */
for (n = 0; n < set_buf.maxis && !axmatch; n++) {
if (axtype == set_buf.ax[n].axtype) { /* matching axis type */
/* does rest match ? */
if (skysys == set_buf.ax[n].skysys && prosys == set_buf.ax[n].prosys && velsys == set_buf.ax[n].velsys) {
axmatch = n + 1; /* matching axis number */
}
}
}
return( axmatch ); /* return to caller */
}
static void addsynonym( fchar cname, fchar ctype )
/*
* addsynonym adds a synonym to syn_buf.
*/
{
fint n; /* loop counter */
for (n = 0; n < cname.l && isprint(cname.a[n]); n++);
if (n != cname.l || cname.a[0] == ' ') return;
/* search loop */
for (n = 0; n < nsyns && fcomp( syn_buf[n].cname, cname ); n++);
if (n == nsyns) { /* increase storage space */
if (nsyns < MAXSYNONYMS) { /* buffer not yet full */
nsyns += 1; /* increase synonym counter */
}
}
if (n < nsyns) { /* can we copy */
fcopy( syn_buf[n].cname, cname ); /* copy synonym */
fcopy( syn_buf[n].ctype, ctype ); /* copy ctype */
}
}
static void getsynonym( fchar cname )
/*
* getsynonym searches the synonym buffer for a synonym for cname.
* If found, it returns the axis name (in cname).
*/
{
fint n; /* loop counter */
/* search loop */
for (n = 0; n < nsyns && fcomp( syn_buf[n].cname, cname ); n++);
if (n < nsyns) { /* found */
fcopy( cname, syn_buf[n].ctype ); /* return synonym */
}
}
static int getdsccode( char *rec )
/*
* getdsccode returns a code which determines what should be done with it.
* 0 means unknown descriptor.
* -1 not a FITS descriptor.
*/
{
int r = 0; /* return value */
if (*rec == ' ') return( -1 ); /* simple solution */
/* search loop */
while ((r < DSCLISTLEN) && strncmp( rec, dsc_buf[r].dsc, strlen( dsc_buf[r].dsc ) ) ) r++;
if (r < DSCLISTLEN) { /* descriptor in list */
if (dsc_buf[r].rep != NULL) {
int k = 0;
while (k < MAXFTSDSCLEN && dsc_buf[r].rep[k]) {
rec[k] = dsc_buf[r].rep[k]; k++;
}
while (k < MAXFTSDSCLEN) rec[k++] = ' ';
}
r = dsc_buf[r].code; /* return code */
} else if (rec[MAXFTSDSCLEN] == '=') {
r = 0; /* record contains descriptor */
} else { /* no descriptor */
r = -1; /* return code */
}
return( r ); /* return to caller */
}
static void getintape( void )
/*
* getintape prompts the user for the name of the input tape device.
* Next it is checked whether the tape device exists by opening it with
* mtopen. If it could be opened getintape is satisfied and returns to
* the caller, otherwise the keyword is cancelled and the user is
* prompted again for the name of the input tape device.
*/
{
fint elev = 4; /* error level (fatal) */
mtid = mtopen_c( INTAPE_DEV ); /* open tape device */
if (mtid < 0) {
error_c( &elev, INTAPE_ERR1 ); /* unknown tape device */
}
(void) mtname_c( &mtid, fts_buf.dev ); /* get name of tape device */
}
static void getinfile( void )
/*
* getinfile prompts the user for the name of the inputfile.
*/
{
int first;
if ( fid == NULL ) {
first = 1;
} else {
first = 0;
(void ) fclose( fid );
}
do {
fint n;
fint dlev = 1;
if ( first ) {
n = usertext_c( fts_buf.dev, &dlev, FITSFILE_KEY, FITSFILE_MES1 );
} else {
n = usertext_c( fts_buf.dev, &dlev, FITSFILE_KEY, FITSFILE_MES2 );
}
if ( n ) {
char name[FILENAME_MAX+1];
strncpy( name, fts_buf.dev.a, n );
name[n] = 0;
fid = fopen( name, "rb" );
if ( fid == NULL ) {
reject_c( FITSFILE_KEY, tofchar( "File not found!" ) );
} else {
nfile = 1;
break;
}
} else {
nfile = 0;
break;
}
} while ( 1 );
cancel_c( FITSFILE_KEY );
}
static void getinfiles( void )
/*
* getinfiles prompts the user for the file numbers on tape. The default
* action is that rfits quits. Negative file numbers mean previous files,
* positive numbers next files and zero the current file. The program is
* setup so that after each load the tape will be positioned after a
* tapemark, so the next read gives us the first header block.
* Only nfile (= number of files) and files[] (= array containing file
* numbers) are exported. The file numbers are sorted in ascending order.
* The keyword is cancelled each time for the repeat loop in the main program.
*/
{
fchar key = INFILES_KEY; /* keyword */
fchar mes = INFILES_MES; /* message */
fint dlev = 1; /* default level */
fint maxfiles = MAXFILES; /* maximum number of files */
if ( fid != NULL ) return; /* ready already */
nfile = userint_c( files, &maxfiles, &dlev, key, mes );
cancel_c( key ); /* cancel keyword */
if (nfile) { /* do not quit yet */
sortia_c( files, &nfile ); /* sort in ascending order */
}
}
static void getoutset( void )
/*
* getoutset prompts the user for the set name. If the set does exist,
* the significant descriptor items of the set are read and stored in
* a structure.
*/
{
fint dlev = 0; /* default level (no default) */
fint gerror = 0; /* for GDS errors */
do {
fchar key = OUTSET_KEY;
fchar mes = OUTSET_MES;
(void) usertext_c( set_buf.set, &dlev, key, mes );
cancel_c( key ); /* cancel keyword */
set_buf.exist = tobool( gds_exist_c( set_buf.set, &gerror ) );
gerror = 0;
if ( set_buf.exist ) {
bool okay = TRUE;
fchar key = OKAY_KEY;
fchar mes = OKAY_MES;
fint dlev = 1;
fint nlog = 1;
(void) userlog_c( &okay, &nlog, &dlev, key, mes );
cancel_c( key );
if ( tobool( okay ) ) break;
} else break;
} while (1);
setfblank_c( &set_buf.blank ); /* set blank value */
if (set_buf.exist) { /* set already exists */
fint level = 0; /* top level */
fint n = 0; /* loop counter */
fint pmask, smask; /* bit masks */
gdsd_rint_c( set_buf.set, tofchar( "NAXIS" ), &level, &set_buf.naxis, &gerror );
if (gerror) gerror = 0;
gdsd_rdble_c( set_buf.set, tofchar( "EPOCH" ), &level, &set_buf.epoch, &gerror );
if (gerror) set_buf.epoch = 0.0; gerror = 0;
gdsd_rdble_c( set_buf.set, tofchar( "FREQ0" ), &level, &set_buf.freq0, &gerror );
if (gerror) set_buf.freq0 = 0.0; gerror = 0;
gdsd_rchar_c( set_buf.set, tofchar( "INSTRUME" ), &level, set_buf.instrume, &gerror );
if (gerror) gerror = 0;
do { /* loop to get axes info */
char ctypeb[MAXFTSNAMLEN], cunitb[MAXFTSNAMLEN];
char dtypeb[MAXFTSNAMLEN], dunitb[MAXFTSNAMLEN];
double cdelt, crota, crpix, crval;
double ddelt, drota, drpix, drval;
fchar ctype, cunit;
fchar dtype, dunit;
fint naxis;
fint skysys, prosys, velsys;
pmask = smask = 0;
ctype = finit( ctypeb, MAXFTSNAMLEN );
cunit = finit( cunitb, MAXFTSNAMLEN );
dtype = finit( dtypeb, MAXFTSNAMLEN );
dunit = finit( dunitb, MAXFTSNAMLEN );
gdsd_rint_c( set_buf.set, descr( "NAXIS", n + 1 ), &level, &naxis, &gerror );
if (gerror) naxis = gerror = 0;
gdsd_rdble_c( set_buf.set, descr( "CDELT", n + 1 ), &level, &cdelt, &gerror );
if (gerror) gerror = 0; else pmask += MDELT;
gdsd_rdble_c( set_buf.set, descr( "CROTA", n + 1 ), &level, &crota, &gerror );
if (gerror) gerror = 0; else pmask += MROTA;
gdsd_rdble_c( set_buf.set, descr( "CRPIX", n + 1 ), &level, &crpix, &gerror );
if (gerror) gerror = 0; else pmask += MRPIX;
gdsd_rdble_c( set_buf.set, descr( "CRVAL", n + 1 ), &level, &crval, &gerror );
if (gerror) gerror = 0; else pmask += MRVAL;
gdsd_rchar_c( set_buf.set, descr( "CTYPE", n + 1 ), &level, ctype, &gerror );
if (gerror) gerror = 0; else pmask += MTYPE;
gdsd_rchar_c( set_buf.set, descr( "CUNIT", n + 1 ), &level, cunit, &gerror );
if (gerror) gerror = 0; else pmask += MUNIT;
gdsd_rdble_c( set_buf.set, descr( "DDELT", n + 1 ), &level, &ddelt, &gerror );
if (gerror) gerror = 0; else smask += MDELT;
gdsd_rdble_c( set_buf.set, descr( "DROTA", n + 1 ), &level, &drota, &gerror );
if (gerror) gerror = 0; else smask += MROTA;
gdsd_rdble_c( set_buf.set, descr( "DRPIX", n + 1 ), &level, &drpix, &gerror );
if (gerror) gerror = 0; else smask += MRPIX;
gdsd_rdble_c( set_buf.set, descr( "DRVAL", n + 1 ), &level, &drval, &gerror );
if (gerror) gerror = 0; else smask += MRVAL;
gdsd_rchar_c( set_buf.set, descr( "DTYPE", n + 1 ), &level, dtype, &gerror );
if (gerror) gerror = 0; else smask += MTYPE;
gdsd_rchar_c( set_buf.set, descr( "DUNIT", n + 1 ), &level, dunit, &gerror );
if (gerror) gerror = 0; else smask += MUNIT;
if (smask || pmask) {
set_buf.ax[n].naxis = naxis;
set_buf.ax[n].pmask = pmask;
set_buf.ax[n].cdelt = cdelt;
set_buf.ax[n].crota = crota;
set_buf.ax[n].crpix = crpix;
set_buf.ax[n].crval = crval;
fcopy( set_buf.ax[n].ctype, ctype );
fcopy( set_buf.ax[n].cunit, cunit );
set_buf.ax[n].smask = smask;
set_buf.ax[n].ddelt = ddelt;
set_buf.ax[n].drota = drota;
set_buf.ax[n].drpix = drpix;
set_buf.ax[n].drval = drval;
fcopy( set_buf.ax[n].dtype, dtype );
fcopy( set_buf.ax[n].dunit, dunit );
set_buf.ax[n].low = 1 - nint( crpix );
set_buf.ax[n].upp = naxis - nint( crpix );
set_buf.ax[n].min = set_buf.ax[n].low;
set_buf.ax[n].max = set_buf.ax[n].upp;
set_buf.ax[n].axtype = axtype_c( ctype, cunit, dunit, &skysys, &prosys, &velsys );
set_buf.ax[n].skysys = skysys; /* sky system code */
set_buf.ax[n].prosys = prosys; /* projection system code */
set_buf.ax[n].velsys = velsys; /* velocity system code */
n += 1; /* increase number of axes */
}
} while (smask || pmask); /* until no information found */
set_buf.maxis = n; /* total number of axis */
}
}
static void position( fint nextfile )
/*
* position positions the tape at the beginning of the file
* denoted by nextfile. Counting is done relative to the initial
* position, which is 0. Negative numbers indicate previous files,
* positive numbers files further on. The static variable currentfile
* holds the current file number. This is the only exported variable.
* It is assumed that the tape is positioned at the first record of a
* file, that is, just after a tapemark or BOT.
*/
{
fint elev = 4; /* error level (FATAL) */
if ( fid != NULL ) return; /* quick exit */
if (nextfile > currentfile) { /* do a forward skip */
fint fsf = nextfile - currentfile; /* number of files to skip */