← Back to branch summary

~ubuntu-branches/ubuntu/wily/eso-midas/wily-proposed

« back to all changes in this revision

Viewing changes to prim/dio/libsrc/fitsinf.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/dio/libsrc/fitsinf.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
.COPYRIGHT  (c)  1995-2009   European Southern Observatory
 
30
.IDENT      fitsinf.c
 
31
.LAUGUAGE   C
 
32
.AUTHOR     P.Grosbol   ESO/IPG
 
33
.KEYWORDS   MIDAS frames, information
 
34
.COMMENT    get informations from MIDAS frames
 
35
.VERSION    1.0  1991-Mar-23 : Creation,   PJG
 
36
.VERSION    1.1  1991-May-15 : Check for SCFGET return,   PJG
 
37
.VERSION    1.2  1991-Jun-03 : Always set bfac and bzero,   PJG
 
38
.VERSION    1.3  1991-Sep-24 : Get NULL's for tables,   PJG
 
39
.VERSION    1.4  1992-Aug-07 : Correct format for C*n table columns, PJG
 
40
.VERSION    1.5  1993-Apr-12 : Correct MIN/MAX computation, PJG
 
41
.VERSION    1.6  1993-Sep-17 : Correct fitstbl because of TCBGET
 
42
.VERSION    1.7  1993-Oct-26 : Update to new TC + prototypes, PJG
 
43
.VERSION    1.8  1995-Jan-23 : Do not update LHCUTS, correct check, PJG
 
44
.VERSION    1.9  1996-Jun-24 : Exchange MAXLONG with 2147483645.0, PJG
 
45
.VERSION    2.0  1998-Mar-09 : use existing BSCALE, BZERO if FITS file already
 
46
.VERSION    2.1  1998-Apr-22 : Use I2 format of I1 table columns, PJG
 
47
 
 
48
 091229         last modif
 
49
---------------------------------------------------------------------*/
 
50
 
 
51
#include   <stdio.h>
 
52
 
 
53
#include   <osparms.h>
 
54
#include   <computer.h>
 
55
#include   <fitsdef.h>
 
56
#include   <fitskwt.h>
 
57
#include   <fileexts.h>
 
58
/*
 
59
#include   <midas_def.h>
 
60
*/
 
61
#include   <tblsys.h>
 
62
 
 
63
#define    BSIZE       10240         /* size of internal buffer      */
 
64
#define    MXLB           81         /* max. char. in line buffer    */
 
65
 
 
66
static     int         fpexc;        /* fp- exception flag           */
 
67
 
 
68
static     SDEF         sdef;        /* Image scaling information    */
 
69
static     TXDEF       txdef;        /* table definitions            */
 
70
static     FDEF    fdef[MXF];        /* information on columns       */
 
71
/*
 
72
 
 
73
*/
 
74
 
 
75
SDEF *fitsbdf(mfd,mff,mfn,outflg)
 
76
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
77
.PURPOSE   compute scaling information MIDAS image
 
78
.COMMENT   if mfd<0 the pointer is returned directly!    
 
79
.RETURN    return pointer to scaling structure, on error NULL pointer
 
80
---------------------------------------------------------------------*/
 
81
int        mfd;                /* IN: MIDAS file number              */
 
82
int        mff;                /* IN: MIDAS file data format         */
 
83
char       *mfn;               /* IN: name of MIDAS file             */
 
84
int        *outflg;            /* OUT: 99, if BSCALE, BZERO are used
 
85
                                       0, if LHCUTS was set or computed */
 
86
{
 
87
int            i, n, ns, nb, nv;
 
88
int            mmm, ioff, null, naxis;
 
89
int            nax[MXDIM], unit[4];
 
90
void           fpeh();
 
91
float          *f, *fb, cuts[2];
 
92
double         lcut, hcut;
 
93
 
 
94
 
 
95
*outflg = 0;
 
96
if (mfd<0) return &sdef;
 
97
 
 
98
SCDRDI(mfd,"NAXIS",1,1,&nv,&naxis,unit,&null);
 
99
sdef.dsize = (naxis) ? 1 : 0;
 
100
SCDRDI(mfd,"NPIX",1,MXDIM,&nv,nax,unit,&null);
 
101
for (i=0; i<naxis; i++) sdef.dsize *= nax[i];
 
102
 
 
103
cuts[0] = cuts[1] = 0.0;
 
104
SCDRDR(mfd,"LHCUTS",3,2,&nv,cuts,unit,&null);
 
105
lcut = cuts[0]; hcut = cuts[1];
 
106
 
 
107
if (nv!=2 || hcut<=lcut) 
 
108
   {                            /* compute max/min cuts of data */
 
109
   struct FCT_STRUCT *fctpntr;
 
110
   int   istat;
 
111
   float  bscale, bzero;
 
112
 
 
113
   fctpntr = FCT.ENTRIES + mfd;
 
114
   if (fctpntr->CR_FLAG == 1)   /* we work on a created FITS file */
 
115
      {                         /* take cuts as they are */
 
116
      sdef.dmin = sdef.dmax = 0.0;
 
117
      sdef.bfac = 1.0;
 
118
      sdef.boff = 0.0;
 
119
      return &sdef;
 
120
      }
 
121
 
 
122
   if (fctpntr->FILTYP > 0) 
 
123
      {                         /* take FITS header as it is */
 
124
      istat = SCDRDR(mfd,"BSCALE",1,1,&i,&bscale,unit,&null);
 
125
      if (istat == ERR_NORMAL) 
 
126
         istat = SCDRDR(mfd,"BZERO",1,1,&i,&bzero,unit,&null);
 
127
      if (istat == ERR_NORMAL) 
 
128
         {
 
129
         sdef.bfac = bscale;
 
130
         sdef.boff = bzero;
 
131
         }
 
132
      else
 
133
         {
 
134
         sdef.bfac = 1.0;
 
135
         sdef.boff = 0.0;
 
136
         }
 
137
 
 
138
      if (nv == 2)              /* in first if-test: low, hi cuts were equal */
 
139
         {
 
140
         sdef.dmin = sdef.dmax = 0.0;
 
141
         }
 
142
      else
 
143
         *outflg = 99;
 
144
 
 
145
      return &sdef;
 
146
      } 
 
147
 
 
148
   fb = (float *)osmmget(BSIZE);        /* get internal buffer */
 
149
   if (!fb) 
 
150
      {
 
151
      SCTPUT("Error: cannot allocate internal buffer");
 
152
      SCFCLO(mfd); return (SDEF *) 0;
 
153
      }
 
154
 
 
155
   osscatch(SIGFPE,fpeh); fpexc = 0;             /* catch fp-exception  */
 
156
 
 
157
   ioff = 1; nb = BSIZE >> 2;                   /* = BSIZE/4 */
 
158
   lcut = hcut = 0.0;
 
159
   if (mff!=D_R4_FORMAT)
 
160
      {                         /* we need data in floating point format */
 
161
      SCFCLO(mfd);
 
162
      SCFOPN(mfn,D_R4_FORMAT,0,F_IMA_TYPE,&mmm);
 
163
      }
 
164
   else 
 
165
      mmm = mfd;
 
166
 
 
167
   while (ioff<=sdef.dsize) 
 
168
      {                       /* find max/min values in data */
 
169
      f = fb;
 
170
      n = (sdef.dsize-ioff<nb) ? sdef.dsize-ioff+1 : nb;
 
171
      SCFGET(mmm, ioff, n, &ns, (char *) f);
 
172
      if (ioff == 1) lcut = hcut = *f;
 
173
      ioff += ns;
 
174
      if (ns<n) ioff = sdef.dsize + 1;    /* to force exit of loop  */
 
175
      while (ns--) 
 
176
         {
 
177
         fpexc = 0;
 
178
         if (!isNULLFLOAT(*f))
 
179
            {
 
180
            if (!fpexc)
 
181
               {
 
182
               if (*f<lcut) lcut = *f;
 
183
               else if (hcut<*f) hcut = *f;
 
184
               }
 
185
            }
 
186
         f++;
 
187
         }
 
188
      }
 
189
 
 
190
   if (mff!=D_R4_FORMAT) 
 
191
      {                                 /* reopen in original data format */
 
192
      SCFCLO(mmm);
 
193
      SCFOPN(mfn,mff,0,F_IMA_TYPE,&mfd);
 
194
      }
 
195
   osmmfree((char *) fb);
 
196
   }
 
197
 
 
198
sdef.dmax = hcut;
 
199
sdef.dmin = lcut;
 
200
sdef.bfac = 0.5 * (hcut-lcut)/2147483645.0;
 
201
if (sdef.bfac==0.0) sdef.bfac = 1.0;
 
202
sdef.boff = 0.5 * (hcut + lcut);
 
203
 
 
204
return &sdef;
 
205
}
 
206
/*
 
207
 
 
208
*/
 
209
 
 
210
TXDEF *fitstbl(mfd,ffmt,cut)
 
211
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
212
.PURPOSE   get table information
 
213
.COMMENT   if mfd<0 the pointer is returned directly!    
 
214
.RETURN    return pointer to Table structure
 
215
---------------------------------------------------------------------*/
 
216
int        mfd;                /* IN: MIDAS file number              */
 
217
int        ffmt;               /* IN: file format  B/O/X             */
 
218
int        cut;                /* IN: cut flag                       */
 
219
{
 
220
  char          *pc, fc;
 
221
  short          snull;
 
222
  int            ns, i, dtype;
 
223
  int            len, fr, fw, fm;
 
224
  int            bytes, items, null;
 
225
  FDEF           *fd;
 
226
 
 
227
 
 
228
 
 
229
  if (mfd<0) return &txdef;
 
230
 
 
231
  TCIGET(mfd,&(txdef.tfields),&(txdef.nrow),&ns,&ns,&ns);
 
232
  if (txdef.tfields > MXF)
 
233
     {
 
234
     char tbuf[80];
 
235
     (void) sprintf(tbuf,"no. of columns (%d) > max. supported columns (%d)",
 
236
                    txdef.tfields,MXF);
 
237
     SCETER(66,tbuf);
 
238
     return 0;                  /* that should always break it... */
 
239
     }
 
240
 
 
241
  txdef.mxrow = 0; txdef.mxcol = 0; txdef.col = fdef;
 
242
  fd = fdef;
 
243
 
 
244
  for (i=1; i<=txdef.tfields; i++, fd++) {   /* get info. for columns */
 
245
    TCFGET(mfd,i,fd->tdisp,&len,&dtype);
 
246
    TCBGET(mfd,i,&dtype,&items,&bytes);
 
247
    pc = fd->tdisp;
 
248
    while (*pc && *pc!=' ') pc++;
 
249
    *pc = '\0';
 
250
    if (ffmt == 'B') fd->trepn =  (dtype!=D_C_FORMAT) ? 1 : bytes;
 
251
    else fd->trepn =  (dtype!=D_C_FORMAT) ? items : bytes;  
 
252
    strcpy(fd->tform,fd->tdisp);
 
253
    dcffmt(fd->tform,&fr,&fc,&fw,&fm);
 
254
    TCLGET(mfd,i,fd->ttype);
 
255
    TCUGET(mfd,i,fd->tunit);
 
256
    fd->tnnul = 0;
 
257
    switch (dtype) {
 
258
       case D_R4_FORMAT :
 
259
            fd->tdfmt = 'E';
 
260
            if (ffmt!='B') fd->twdth = 4*items;
 
261
            else if (cut=='C') {
 
262
              if (fc=='I') fd->tdfmt = 'I';
 
263
              fd->twdth = fw;
 
264
            }
 
265
            else { strcpy(fd->tform,"E15.8"); fd->twdth = 15; }
 
266
            break;
 
267
       case D_R8_FORMAT :
 
268
            fd->tdfmt = 'D';
 
269
            if (ffmt!='B') fd->twdth = 8*items;
 
270
            else if (cut=='C') fd->twdth = fw;
 
271
            else {
 
272
              strcpy(fd->tform,"D24.16"); fd->twdth = 24;
 
273
            }
 
274
            break;
 
275
       case D_C_FORMAT :
 
276
            fd->tdfmt = 'A';
 
277
            if (ffmt!='B') fd->twdth = bytes;
 
278
            else if (cut=='C') fd->twdth = fw;
 
279
            else {
 
280
              sprintf(fd->tform,"A%d",bytes); fd->twdth = bytes;
 
281
            }
 
282
            break;
 
283
       case D_I1_FORMAT :
 
284
            fd->tdfmt = 'S';
 
285
            if (ffmt!='B') fd->twdth = 2*items;
 
286
            else if (cut=='C') fd->twdth = fw;
 
287
            else {
 
288
              strcpy(fd->tform,"I4"); fd->twdth = 4;
 
289
            }
 
290
            TBL_toNULL((TBL_D_I2<<TBL_D_BITS)|1, (char *) &snull);
 
291
            fd->tnnul = snull;
 
292
            break;
 
293
       case D_I2_FORMAT :
 
294
            fd->tdfmt = 'S';
 
295
            if (ffmt!='B') fd->twdth = 2*items;
 
296
            else if (cut=='C') fd->twdth = fw;
 
297
            else {
 
298
              strcpy(fd->tform,"I6"); fd->twdth = 6;
 
299
            }
 
300
            TBL_toNULL((TBL_D_I2<<TBL_D_BITS)|1, (char *) &snull);
 
301
            fd->tnnul = snull;
 
302
            break;
 
303
       case D_I4_FORMAT :
 
304
            fd->tdfmt = 'I';
 
305
            if (ffmt!='B') fd->twdth = 4*items;
 
306
            else if (cut=='C') fd->twdth = fw;
 
307
            else {
 
308
              strcpy(fd->tform,"I11"); fd->twdth = 11;
 
309
            }
 
310
            TBL_toNULL((TBL_D_I4<<TBL_D_BITS)|1, (char *) &null);
 
311
            fd->tnnul = null;
 
312
            break;
 
313
     }
 
314
    txdef.mxrow += fd->twdth;
 
315
    if (txdef.mxcol<fd->twdth) txdef.mxcol = fd->twdth;
 
316
  }
 
317
 
 
318
  return &txdef;
 
319
}
 
320
 
 
321
 
 
322
 
 
323
void fpeh(s)                                 /* fp-exception handler   */
 
324
int  s;
 
325
{
 
326
osscatch(SIGFPE,fpeh);
 
327
fpexc++;
 
328
}