← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/dio/libsrc/fitsrdmNUL.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/fitsrdmNUL.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1994-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
.IDENT      fitsrdm.c
 
30
.LAUGUAGE   C
 
31
.AUTHOR     P.Grosbol   ESO/IPG
 
32
.KEYWORDS   FITS data matrix, decode, read
 
33
.COMMENT    read the prime data matrix of FITS file
 
34
.VERSION    1.0  1988-Dec-10 : Creation,   PJG 
 
35
.VERSION    1.1  1989-Oct-08 : Save Group parm. in table,   PJG 
 
36
.VERSION    1.2  1989-Oct-24 : Check of data to store,   PJG 
 
37
.VERSION    1.3  1990-Feb-04 : Change call-seq. for cv-routine, PJG 
 
38
.VERSION    1.4  1991-Feb-14 : Correct pointer error, PJG
 
39
.VERSION    2.0  1991-Feb-24 : Change call-seq. structures, PJG 
 
40
.VERSION    2.1  1991-Mar-23 : Compute min/max if not given, PJG
 
41
.VERSION    2.2  1993-Sep-23 : Check overflow in cut values, PJG
 
42
.VERSION    2.3  1993-Dec-02 : Update to new SC + prototypes, PJG
 
43
.VERSION    2.4  1994-May-09 : Cast variable for ANSI-C, PJG
 
44
.VERSION    2.5  1994-Jun-28 : Include UI2 as bitpix=-16, PJG
 
45
.VERSION    2.6  1994-Sep-07 : Add explicit cast for unsigned, PJG
 
46
.VERSION    2.7  1997-Jan-08 : Only warning on wrong record size, PJG
 
47
 
 
48
 090706         last modif
 
49
---------------------------------------------------------------------*/
 
50
 
 
51
#include   <stdio.h>            /* computer specific constants    */
 
52
 
 
53
#include   <computer.h>            /* computer specific constants    */
 
54
#include   <fitsfmt.h>             /* general data definitions       */
 
55
#include   <fitsdef.h>             /* basic FITS definitions         */
 
56
#include   <fitsextvals.h>
 
57
#include   <midas_def.h>           /* MIDAS definitions              */
 
58
 
 
59
#define    MXFB              2880  /* max. size of scaling buffer    */
 
60
 
 
61
/*
 
62
 
 
63
*/
 
64
 
 
65
#ifdef __STDC__
 
66
int fitsrdmNUL(int mfd , BFDEF * bfdef , int size , int mfdt , char fmt,
 
67
            int Midas_flag)
 
68
#else
 
69
int fitsrdmNUL(mfd,bfdef,size,mfdt,fmt,Midas_flag)
 
70
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
71
.PURPOSE        Read prime data matrix in FITS
 
72
.ALGORITHM      this version does the NULL check of data
 
73
.RETURN         error code - 0: OK, -1: error
 
74
---------------------------------------------------------------------*/
 
75
int       mfd;                 /* IN: MIDAS file descriptor          */
 
76
BFDEF     *bfdef;              /* IN: pointer to FITS header parm's  */
 
77
int       size;                /* IN: size of data matrix (bytes)    */
 
78
int       mfdt;                /* IN: MIDAS file desc. Groups table  */
 
79
char      fmt;                 /* IN: data format Original/Float     */
 
80
int       Midas_flag;          /* IN: = 0, `usual' mode
 
81
                                      = 1, `internal FITS access'  */
 
82
 
 
83
#endif
 
84
{
 
85
register int   bf, k, ireg;
 
86
unsigned char  uc, *pb;
 
87
short          s, *ps;
 
88
unsigned short u, *pu;
 
89
int            fpf, dno, gno, pno, npix, n, fmm;
 
90
int            dn, pn, ndata, dfmt, nb, pcnt, dsf, unit[4];
 
91
int            felem, i, *pi, mnvi, mxvi;
 
92
float          *pf, cuts[2];
 
93
register double      fac, zero, *pd, mnvd, mxvd;
 
94
double         d, dreg;
 
95
PDEF           *pp;
 
96
union { 
 
97
        float    f[2*MXFB];
 
98
        double     d[MXFB];
 
99
      } buf;
 
100
union { 
 
101
        unsigned char   *b;
 
102
        char            *c;
 
103
        short           *s;
 
104
        unsigned short  *u;
 
105
        int             *i;
 
106
        float           *f;
 
107
        double          *d;
 
108
      } p;
 
109
 
 
110
 
 
111
int    ifac;
 
112
 
 
113
 
 
114
mnvi = mxvi = 0;
 
115
mnvd = mxvd = 0.0;
 
116
 
 
117
pno = bfdef->pcount; pcnt = 0;             /* initiate counters     */
 
118
pp = bfdef->parm;
 
119
npix = gno = 0;
 
120
bf = bfdef->bflag;
 
121
dfmt = bfdef->bitpix;
 
122
nb = (dfmt<0) ? -dfmt/8 : dfmt/8;
 
123
ndata = size/(nb*bfdef->gcount) - pno; dno = ndata;
 
124
fac = bfdef->bscale; zero = bfdef->bzero;
 
125
if ((fac > 0.999999) && (fac < 1.00001))
 
126
   ifac = 1;                            /* no need to multiply with 1.0 */
 
127
else
 
128
   ifac = 0;
 
129
felem = 1;
 
130
fmm = ((bfdef->mflag != 3) && (0<size));
 
131
fpf = (bfdef->sflag || dfmt!=-32) ;
 
132
dsf = ((bfdef->sflag && dfmt!=-64) || fmt=='F') ? -32 : dfmt;
 
133
 
 
134
while (0<size) 
 
135
   {                 /* read all data in prime matrix  */
 
136
   if ((n=dread(&p.c,FITSLR)) != FITSLR)        /* read next data record  */
 
137
      {                 
 
138
      if (size<=n) 
 
139
         SCTPUT("Warning: incomplete FITS record read!");
 
140
      else 
 
141
         {
 
142
         int   unit;
 
143
         char  tbuf[80];
 
144
 
 
145
         SCTPUT("Error: unexpected EOF");
 
146
         switch (dfmt)
 
147
            {
 
148
            case   8 : break;
 
149
            case -16 :
 
150
            case  16 : size /= 2; break;
 
151
            case  32 : 
 
152
            case -32 : size /= 4; break;
 
153
            case -64 : size /= 8; break;
 
154
            }
 
155
         (void) sprintf(tbuf,"%d data values still missing",size);
 
156
         SCTPUT(tbuf);
 
157
         if (0<=mfd) SCFCLO(mfd); mfd = -1;
 
158
         (void) SCKWRI("OUTPUTI",&size,16,1,&unit);     /* save value */
 
159
         return NOFITS;
 
160
         }
 
161
      }
 
162
 
 
163
   if (size<=n) 
 
164
      {
 
165
      n = size;                 /* last loop */
 
166
      size = 0;
 
167
      }
 
168
   else
 
169
      size -= n;                /* decrement remaining bytes      */
 
170
 
 
171
   switch (dfmt) 
 
172
      {                /* convert to local data format   */
 
173
      case   8 : npix = n; break;
 
174
      case -16 :
 
175
      case  16 : npix = n/2; 
 
176
                 if (!same_comp_i2) cvi2(p.s,npix,0); break;
 
177
      case  32 : npix = n/4; 
 
178
                 if (!same_comp_i4) cvi4(p.i,npix,0); break;
 
179
      case -32 : npix = n/4; cvr4(p.f,npix,0); break;
 
180
      case -64 : npix = n/8; cvr8(p.d,npix,0); break;
 
181
      }
 
182
 
 
183
   do 
 
184
      {                         /* scale all values if needed     */
 
185
      if (0<pno) 
 
186
         {                      /* decode groups parameters       */
 
187
         pn = (pno<npix) ? pno : npix;
 
188
         pno -= pn; npix -= pn;
 
189
         while (pn--) 
 
190
            {
 
191
            switch (dfmt) 
 
192
               {
 
193
               case   8 : d = pp->pscal * (*(p.b++)) + pp->pzero;
 
194
                          break;
 
195
               case -16 :
 
196
               case  16 : d = pp->pscal * (*(p.s++)) + pp->pzero;
 
197
                          break;
 
198
               case  32 : d = pp->pscal * (*(p.i++)) + pp->pzero;
 
199
                          break;
 
200
               case -32 : d = pp->pscal * (*(p.f++)) + pp->pzero;
 
201
                          break;
 
202
               case -64 : d = pp->pscal * (*(p.d++)) + pp->pzero;
 
203
                          break;
 
204
               }
 
205
            pcnt++; pp++;
 
206
            if (0<=mfdt) TCEWRD(mfdt,gno+1,pcnt,&d);
 
207
            }
 
208
         }
 
209
 
 
210
      if (!pno && 0<dno && npix) 
 
211
         {                              /* decode data values            */
 
212
         dn = (dno<npix) ? dno : npix;
 
213
         dno -= dn; npix -= dn;
 
214
 
 
215
         if (bfdef->sflag) 
 
216
            {                           /* we scale the data */
 
217
            switch (dfmt) 
 
218
               {
 
219
               case   8 : 
 
220
                  pf = buf.f; uc = bfdef->blank;
 
221
                  if (bf)
 
222
                     {
 
223
                     if (ifac == 0)
 
224
                        {
 
225
                        for (k=0; k<dn; k++, pf++, p.b++) 
 
226
                           {
 
227
                           if (*p.b==uc) 
 
228
                              toNULLFLOAT(*pf);
 
229
                           else
 
230
                              *pf = fac*((int)(*p.b)) + zero;
 
231
                           }
 
232
                        }
 
233
                     else               /* fac = 1.0 */
 
234
                        {
 
235
                        for (k=0; k<dn; k++, pf++, p.b++) 
 
236
                           {
 
237
                           if (*p.b==uc) 
 
238
                              toNULLFLOAT(*pf);
 
239
                           else
 
240
                              *pf = (int)(*p.b) + zero;
 
241
                           }
 
242
                        }
 
243
                     }
 
244
                  else
 
245
                     {
 
246
                     if (ifac == 0)
 
247
                        {
 
248
                        for (k=0; k<dn; k++)
 
249
                            *pf++ = fac*((int)(*p.b++)) + zero;
 
250
                        }
 
251
                     else               /* fac = 1.0 */
 
252
                        {
 
253
                        for (k=0; k<dn; k++) *pf++ = (int)(*p.b++) + zero;
 
254
                        }
 
255
                     }
 
256
                  break;
 
257
               case -16 :
 
258
               case  16 :
 
259
                  pf = buf.f; s = bfdef->blank;
 
260
                  if (bf)
 
261
                     {
 
262
                     if (ifac == 0)
 
263
                        {
 
264
                        for (k=0; k<dn; k++, pf++, p.s++) 
 
265
                           {
 
266
                           if (*p.s==s) 
 
267
                              toNULLFLOAT(*pf);
 
268
                           else
 
269
                              *pf = fac*(*p.s) + zero;
 
270
                           }
 
271
                        }
 
272
                     else
 
273
                        {
 
274
                        for (k=0; k<dn; k++, pf++, p.s++) 
 
275
                           {
 
276
                           if (*p.s==s) 
 
277
                              toNULLFLOAT(*pf);
 
278
                           else
 
279
                              *pf = *p.s + zero;
 
280
                           }
 
281
                        }
 
282
                     }
 
283
                  else
 
284
                     {
 
285
                     if (ifac == 0)
 
286
                        {
 
287
                        for (k=0; k<dn; k++)
 
288
                           *pf++ = fac*(*p.s++) + zero;
 
289
                        }
 
290
                     else
 
291
                        {
 
292
                        for (k=0; k<dn; k++) *pf++ = (*p.s++) + zero;
 
293
                        }
 
294
                     }
 
295
 
 
296
                  break;
 
297
               case  32 : 
 
298
                  pf = buf.f; i = bfdef->blank;
 
299
                  if (bf)
 
300
                     {
 
301
                     if (ifac == 0)
 
302
                        {
 
303
                        for (k=0; k<dn; k++, pf++, p.i++) 
 
304
                           {
 
305
                           if (*p.i==i) 
 
306
                              toNULLFLOAT(*pf);
 
307
                           else
 
308
                              *pf = fac*(*p.i) + zero;
 
309
                           }
 
310
                        }
 
311
                     else
 
312
                        {
 
313
                        for (k=0; k<dn; k++, pf++, p.i++)
 
314
                           {
 
315
                           if (*p.i==i)
 
316
                              toNULLFLOAT(*pf);
 
317
                           else
 
318
                              *pf = (*p.i) + zero;
 
319
                           }
 
320
                        }
 
321
                     }
 
322
                  else
 
323
                     {
 
324
                     if (ifac == 0)
 
325
                        {
 
326
                        for (k=0; k<dn; k++) *pf++ = fac*(*p.i++) + zero;
 
327
                        }
 
328
                     else
 
329
                        {
 
330
                        for (k=0; k<dn; k++) *pf++ = (*p.i++) + zero;
 
331
                        }
 
332
                     }
 
333
                  break;
 
334
               case -32 : 
 
335
                  pf = buf.f;
 
336
                  if (ifac == 0)
 
337
                     {
 
338
                     for (k=0; k<dn; k++) *pf++ = fac * (*p.f++) + zero;
 
339
                     }
 
340
                  else
 
341
                     {
 
342
                     for (k=0; k<dn; k++) *pf++ = (*p.f++) + zero;
 
343
                     }
 
344
                  break;
 
345
               case -64 :
 
346
                  pd = buf.d;
 
347
                  if (ifac == 0)
 
348
                     {
 
349
                     for (k=0; k<dn; k++) *pd++ = fac * (*p.d++) + zero;
 
350
                     }
 
351
                  else
 
352
                     {
 
353
                     for (k=0; k<dn; k++) *pd++ = (*p.d++) + zero;
 
354
                     }
 
355
                  break;
 
356
               }                /* end switch */
 
357
 
 
358
            if (Midas_flag != 0)
 
359
               {
 
360
               if (dfmt!=-64)             /* store data in MIDAS file */
 
361
                  MyPut(-32,felem,dn,(char *) buf.f);
 
362
               else  
 
363
                  MyPut(dfmt,felem,dn,(char *) buf.d);
 
364
               }
 
365
            else
 
366
               {
 
367
               if (dfmt!=-64)             /* store data in MIDAS file */
 
368
                  SCFPUT(mfd,felem,dn,(char *) buf.f);
 
369
               else
 
370
                  SCFPUT(mfd,felem,dn,(char *) buf.d);
 
371
               }
 
372
            }
 
373
 
 
374
         else if (fmt=='F') 
 
375
            {                           /* we convert data to float */
 
376
            switch (dfmt) 
 
377
               {
 
378
               case   8 : 
 
379
                  pf = buf.f; uc = bfdef->blank;
 
380
                  if (bf)
 
381
                     {
 
382
                     for (k=0; k<dn; k++, pf++, p.b++)
 
383
                        {
 
384
                        if (*p.b==uc)
 
385
                           toNULLFLOAT(*pf);
 
386
                        else
 
387
                           *pf = *p.b;
 
388
                        }
 
389
                     }
 
390
                  else
 
391
                     {
 
392
                     for (k=0; k<dn; k++) *pf++ = *p.b++;
 
393
                     }
 
394
                  break;
 
395
               case -16 :
 
396
               case  16 :
 
397
                  pf = buf.f; s = bfdef->blank;
 
398
                  if (bf)
 
399
                     {
 
400
                     for (k=0; k<dn; k++, pf++, p.s++)
 
401
                        {
 
402
                        if (*p.s==s)
 
403
                           toNULLFLOAT(*pf);
 
404
                        else
 
405
                           *pf = *p.s;
 
406
                        }
 
407
                     }
 
408
                  else
 
409
                     {
 
410
                     for (k=0; k<dn; k++) *pf++ = *p.s++;
 
411
                     }
 
412
 
 
413
                  break;
 
414
               case  32 :
 
415
                  pf = buf.f; i = bfdef->blank;
 
416
                  if (bf)
 
417
                     {
 
418
                     for (k=0; k<dn; k++, pf++, p.i++)
 
419
                        {
 
420
                        if (*p.i==i)
 
421
                           toNULLFLOAT(*pf);
 
422
                        else
 
423
                           *pf = *p.i;
 
424
                        }
 
425
                     }
 
426
                  else
 
427
                     {
 
428
                     for (k=0; k<dn; k++) *pf++ = *p.i++;
 
429
                     }
 
430
                  break;
 
431
               case -32 : 
 
432
                  if (Midas_flag != 0)
 
433
                     MyPut(dfmt,felem,dn,(char *) p.f);
 
434
                  else
 
435
                     SCFPUT(mfd,felem,dn,(char *) p.f);
 
436
                  p.f += dn;
 
437
                  break;
 
438
               case -64 :
 
439
                  pf = buf.f;
 
440
                  for (k=0; k<dn; k++, pf++, p.d++) 
 
441
                     {
 
442
                     if (isNULLDOUBLE(*p.d)) 
 
443
                        toNULLFLOAT(*pf);
 
444
                     else 
 
445
                        *pf = *p.d;
 
446
                     }
 
447
                  break;
 
448
               }                        /* end switch */
 
449
            if (Midas_flag != 0)
 
450
               {
 
451
               if (dfmt!=-32)             /* store data in MIDAS file */
 
452
                  MyPut(-32,felem,dn,(char *) buf.f);
 
453
               }
 
454
            else
 
455
               {
 
456
               if (dfmt!=-32)             /* store data in MIDAS file */
 
457
                  SCFPUT(mfd,felem,dn,(char *) buf.f);
 
458
               }
 
459
            }
 
460
 
 
461
         else 
 
462
            {                           /* just check for NULLs */
 
463
            switch (dfmt)
 
464
               {
 
465
               case   8 : 
 
466
                  if (bf) 
 
467
                     { 
 
468
                     pb = p.b; uc = bfdef->blank;
 
469
                     for (k=0; k<dn; k++, pb++) if (*pb==uc) *pb = 0xFF;
 
470
                     }
 
471
                  if (Midas_flag != 0)
 
472
                     MyPut(dfmt,felem,dn,(char *) p.b);
 
473
                  else
 
474
                     SCFPUT(mfd,felem,dn,(char *) p.b);
 
475
                  p.b += dn;
 
476
                  break;
 
477
               case -16 : 
 
478
                  pu = p.u; u = bfdef->blank; ps = p.s;
 
479
                  if (bf)
 
480
                     {
 
481
                     for (k=0; k<dn; k++, pu++, ps++) 
 
482
                        {
 
483
                        if (*ps==u) 
 
484
                           toNULLSHORT(*pu);
 
485
                        else
 
486
                           *pu = *ps + 32768;
 
487
                        }
 
488
                     }
 
489
                  else
 
490
                     {
 
491
                     for (k=0; k<dn; k++) *pu++ = *ps++ + 32768;
 
492
                     }
 
493
                  if (Midas_flag != 0)
 
494
                     MyPut(dfmt,felem,dn,(char *) p.u);
 
495
                  else
 
496
                     SCFPUT(mfd,felem,dn,(char *) p.u);
 
497
                  p.u += dn;
 
498
                  break;
 
499
               case  16 :
 
500
                  if (bf) 
 
501
                     {
 
502
                     ps = p.s; s = bfdef->blank;
 
503
                     for (k=0; k<dn; k++, ps++)
 
504
                        if (*ps==s) toNULLSHORT(*ps);
 
505
                     }
 
506
                  if (Midas_flag != 0)
 
507
                     MyPut(dfmt,felem,dn,(char *) p.s);
 
508
                  else
 
509
                     SCFPUT(mfd,felem,dn,(char *) p.s);
 
510
                  p.s += dn;
 
511
                  break;
 
512
               case  32 : 
 
513
                  if (bf) 
 
514
                     {
 
515
                     pi = p.i; i = bfdef->blank;
 
516
                     for (k=0; k<dn; k++, pi++)
 
517
                        if (*pi==i) toNULLINT(*pi);
 
518
                     }
 
519
                  if (Midas_flag != 0)
 
520
                     MyPut(dfmt,felem,dn,(char *) p.i);
 
521
                  else
 
522
                     SCFPUT(mfd,felem,dn,(char *) p.i);
 
523
                  p.i += dn;
 
524
                  break;
 
525
               case -32 : 
 
526
                  if (Midas_flag != 0)
 
527
                     MyPut(dfmt,felem,dn,(char *) p.f);
 
528
                  else
 
529
                     SCFPUT(mfd,felem,dn,(char *) p.f);
 
530
                  p.f += dn;
 
531
                  break;
 
532
               case -64 : 
 
533
                  if (Midas_flag != 0)
 
534
                     MyPut(dfmt,felem,dn,(char *) p.d);
 
535
                  else
 
536
                     SCFPUT(mfd,felem,dn,(char *) p.d);
 
537
                  p.d += dn;
 
538
                  break;
 
539
               }                /* end switch */
 
540
            }
 
541
 
 
542
 
 
543
         /* for all 3 cases, find data min/max values if not known 
 
544
            this means another pass over the data */
 
545
 
 
546
         if (fmm) 
 
547
            {    
 
548
            switch (dsf) 
 
549
               {
 
550
               case   8 :
 
551
                  pb = p.b - dn;
 
552
                  if (felem==1) mnvi = mxvi = *pb;
 
553
                  for (k=0; k<dn; k++)
 
554
                     {
 
555
                     ireg = (int) *pb++;
 
556
                     if (ireg != 0xFF) 
 
557
                        {
 
558
                        if (ireg<mnvi) 
 
559
                              mnvi = ireg;
 
560
                        else if (mxvi<ireg) 
 
561
                              mxvi = ireg;
 
562
                        }
 
563
                     }
 
564
                  break;
 
565
               case -16 : 
 
566
                  pu = p.u - dn;
 
567
                  if (felem==1) mnvi = mxvi = *pu;
 
568
                  for (k=0; k<dn; k++)
 
569
                     {
 
570
                     ireg = (int) *pu++;
 
571
                     if (!isNULLSHORT(ireg)) 
 
572
                        {
 
573
                        if (ireg < mnvi)
 
574
                                  mnvi = ireg;
 
575
                        else if (mxvi<ireg)
 
576
                                  mxvi = ireg;
 
577
                        }
 
578
                     }
 
579
                  break;
 
580
               case  16 : 
 
581
                  ps = p.s - dn;
 
582
                  if (felem==1) mnvi = mxvi = *ps;
 
583
                  for (k=0; k<dn; k++)
 
584
                     {
 
585
                     ireg = (int) *ps++;
 
586
                     if (!isNULLSHORT(ireg)) 
 
587
                        {
 
588
                        if (ireg<mnvi) 
 
589
                           mnvi = ireg;
 
590
                        else if (mxvi<ireg) 
 
591
                           mxvi = ireg;
 
592
                        }
 
593
                     }
 
594
                  break;
 
595
               case  32 : 
 
596
                  pi = p.i - dn;
 
597
                  if (felem==1) mnvi = mxvi = *pi;
 
598
                  for (k=0; k<dn; k++)
 
599
                     {
 
600
                     ireg = (int) *pi++;
 
601
                     if (!isNULLINT(ireg)) 
 
602
                        {
 
603
                        if (ireg<mnvi) 
 
604
                           mnvi = ireg;
 
605
                        else if (mxvi<ireg) 
 
606
                           mxvi = ireg;
 
607
                        }
 
608
                     }
 
609
                  break;
 
610
               case -32 : 
 
611
                  pf = (fpf) ? buf.f : p.f - dn;
 
612
                  if (felem==1) mnvd = mxvd = *pf;
 
613
                  for (k=0; k<dn; k++)
 
614
                     {
 
615
                     dreg = (double) *pf++;
 
616
                     if (!isNULLFLOAT(dreg)) 
 
617
                        {
 
618
                        if (dreg<mnvd) 
 
619
                           mnvd = dreg;
 
620
                        else if (mxvd<dreg) 
 
621
                           mxvd = dreg;
 
622
                        }
 
623
                     }
 
624
                  break;
 
625
               case -64 : 
 
626
                  pd = (bfdef->sflag) ? buf.d : p.d - dn;
 
627
                  if (felem==1) mnvd = mxvd = *pd;
 
628
                  for (k=0; k<dn; k++)
 
629
                     {
 
630
                     dreg = *pd++;
 
631
                     if (!isNULLDOUBLE(dreg))
 
632
                        {
 
633
                        if (dreg<mnvd) 
 
634
                           mnvd = dreg;
 
635
                        else if (mxvd<dreg) 
 
636
                           mxvd = dreg;
 
637
                        }
 
638
                     }
 
639
                  break;
 
640
               }                /* end switch */
 
641
            }
 
642
 
 
643
         felem += dn;
 
644
         if (!dno) 
 
645
            {
 
646
            gno++;
 
647
            pno = bfdef->pcount; pcnt = 0;
 
648
            pp = bfdef->parm;
 
649
            dno = ndata;
 
650
            }
 
651
         }
 
652
 
 
653
      } while (npix && gno<bfdef->gcount);
 
654
 
 
655
   }
 
656
 
 
657
if (fmm) 
 
658
   {                         /* save computed min/max values  */
 
659
   if (dsf > -32) 
 
660
      {                         /* handle also UI2 here */
 
661
      cuts[0] = (float)mnvi; cuts[1] = (float)mxvi;
 
662
      }
 
663
   else 
 
664
      {
 
665
      if (MAXFLOAT<mnvd) 
 
666
         mnvd = MAXFLOAT;
 
667
      else if (mnvd<MINFLOAT) 
 
668
         mnvd = MINFLOAT;
 
669
      if (MAXFLOAT<mxvd) 
 
670
         mxvd = MAXFLOAT;
 
671
      else if (mxvd<MINFLOAT) 
 
672
         mxvd = MINFLOAT;
 
673
      cuts[0] = (float)mnvd; cuts[1] = (float)mxvd;
 
674
      }
 
675
    SCDWRR(mfd,"LHCUTS",cuts,3,2,unit);
 
676
 
 
677
   }
 
678
 
 
679
 
 
680
if (Midas_flag == 0)
 
681
   {
 
682
   if (0<=mfd) SCFCLO(mfd);                       /* close data files  */
 
683
   if (0<=mfdt) { TCSINI(mfdt); TCTCLO(mfdt); }
 
684
   mfd = -1; mfdt = -1;
 
685
   }
 
686
 
 
687
return 0;
 
688
}
 
689
/*
 
690
 
 
691
*/
 
692
 
 
693
#ifdef __STDC__
 
694
int fitsrdmUI2NUL(int mfd , BFDEF * bfdef , int size , int mfdt , char fmt,
 
695
            int Midas_flag)
 
696
#else
 
697
int fitsrdmUI2NUL(mfd,bfdef,size,mfdt,fmt,Midas_flag)
 
698
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
699
.PURPOSE       Read prime unsinged short int data matrix in FITS
 
700
.RETURN        error code - 0: OK, -1: error
 
701
---------------------------------------------------------------------*/
 
702
int       mfd;                 /* IN: MIDAS file descriptor          */
 
703
BFDEF     *bfdef;              /* IN: pointer to FITS header parm's  */
 
704
int       size;                /* IN: size of data matrix (bytes)    */
 
705
int       mfdt;                /* IN: MIDAS file desc. Groups table  */
 
706
char      fmt;                 /* IN: data format Original/Float     */
 
707
int       Midas_flag;          /* IN: = 0, `usual' mode
 
708
                                      = 1, `internal FITS access'  */
 
709
 
 
710
#endif
 
711
{
 
712
register int   bf, k, ireg;
 
713
short          s, *ps;
 
714
unsigned short u, *pu;
 
715
int            dno, gno, pno, npix, n, fmm;
 
716
int            dn, pn, ndata, dfmt, pcnt, dsf, unit[4];
 
717
int            ifac, felem, mnvi, mxvi;
 
718
 
 
719
float          *pf, cuts[2];
 
720
 
 
721
double         d;
 
722
register double     dreg, fac, zero, mnvd, mxvd;
 
723
 
 
724
PDEF           *pp;
 
725
union { 
 
726
        float    f[2*MXFB];
 
727
        double     d[MXFB];
 
728
      } buf;
 
729
union { 
 
730
        unsigned char   *b;
 
731
        char            *c;
 
732
        short           *s;
 
733
        unsigned short  *u;
 
734
        int             *i;
 
735
        float           *f;
 
736
        double          *d;
 
737
      } p;
 
738
 
 
739
 
 
740
 
 
741
 
 
742
 
 
743
 
 
744
/* we only come here, if bfdef->bitpix = 16 (or -16) */
 
745
 
 
746
mnvi = mxvi = 0;
 
747
mnvd = mxvd = 0.0;
 
748
pno = bfdef->pcount; pcnt = 0;             /* initiate counters     */
 
749
pp = bfdef->parm;
 
750
gno = 0;
 
751
bf = bfdef->bflag;
 
752
dfmt = bfdef->bitpix;
 
753
/* was: nb = (dfmt<0) ? -dfmt/8 : dfmt/8;
 
754
ndata = size/(nb*bfdef->gcount) - pno; dno = ndata; */
 
755
/* now: since nb = 2 */
 
756
ndata = size/(2*bfdef->gcount) - pno; dno = ndata;
 
757
fac = bfdef->bscale; zero = bfdef->bzero;
 
758
if ((fac > 0.999999) && (fac < 1.00001))
 
759
   ifac = 1;                            /* no need to multiply with 1.0 */
 
760
else
 
761
   ifac = 0;
 
762
felem = 1;
 
763
fmm = ((bfdef->mflag != 3) && (0<size));
 
764
dsf = (bfdef->sflag || fmt=='F') ? -32 : dfmt;
 
765
 
 
766
/*
 
767
printf("UI2: dfmt = %d, bf = %d, fmm = %d, dsf = %d, pno = %d\n",
 
768
dfmt,bf,fmm,dsf,pno);
 
769
*/
 
770
 
 
771
 
 
772
while (0<size) 
 
773
   {                 /* read all data in prime matrix  */
 
774
   if ((n=dread(&p.c,FITSLR)) != FITSLR)        /* read next data record  */
 
775
      {                 
 
776
      if (size<=n) 
 
777
         SCTPUT("Warning: incomplete FITS record read!");
 
778
      else 
 
779
         {
 
780
         int   unit;
 
781
         char  tbuf[80];
 
782
 
 
783
         SCTPUT("Error: unexpected EOF");
 
784
         size /= 2; 
 
785
         (void) sprintf(tbuf,"%d data values still missing",size);
 
786
         SCTPUT(tbuf);
 
787
         if (0<=mfd) SCFCLO(mfd); mfd = -1;
 
788
         (void) SCKWRI("OUTPUTI",&size,16,1,&unit);     /* save value */
 
789
         return NOFITS;
 
790
         }
 
791
      }
 
792
 
 
793
   if (size<n) n = size;
 
794
   size -= n;                     /* decrement remaining bytes      */
 
795
   npix = n/2; 
 
796
   if (!same_comp_i2) cvi2(p.s,npix,0); 
 
797
 
 
798
   do 
 
799
      {                         /* scale all values if needed     */
 
800
      if (0<pno) 
 
801
         {                      /* decode groups parameters       */
 
802
         pn = (pno<npix) ? pno : npix;
 
803
         pno -= pn; npix -= pn;
 
804
         while (pn--) 
 
805
            {
 
806
            d = pp->pscal * (*(p.s++)) + pp->pzero;
 
807
            pcnt++; pp++;
 
808
            if (0<=mfdt) TCEWRD(mfdt,gno+1,pcnt,&d);
 
809
            }
 
810
         }
 
811
 
 
812
      if (!pno && 0<dno && npix) 
 
813
         {                              /* decode data values            */
 
814
         dn = (dno<npix) ? dno : npix;
 
815
         dno -= dn; npix -= dn;
 
816
 
 
817
         /* option = scale the data (unsigned int and int are equal) */
 
818
         /* format for minmax is float */
 
819
 
 
820
         if (bfdef->sflag) 
 
821
            {           
 
822
            pf = buf.f; s = bfdef->blank, ps = p.s;
 
823
            if (fmm)
 
824
               {                        /* dsf = -32 */
 
825
               if (felem==1) mnvd = mxvd = (fac*(*ps)) + zero;
 
826
               if (bf)
 
827
                  { 
 
828
                  for (k=0; k<dn; k++, pf++, ps++)
 
829
                     {
 
830
                     if (*ps == s) 
 
831
                        toNULLFLOAT(*pf);
 
832
                     else
 
833
                        {
 
834
                        *pf = (fac * (*ps)) + zero;
 
835
                        dreg = (double) *pf;
 
836
                        if (dreg<mnvd)
 
837
                           mnvd = dreg;
 
838
                        else if (mxvd<dreg)
 
839
                           mxvd = dreg;
 
840
                        }
 
841
                     }
 
842
                  }
 
843
               else
 
844
                  {
 
845
                  if (ifac == 0)
 
846
                     {
 
847
                     for (k=0; k<dn; k++)
 
848
                        {
 
849
                        *pf = fac*(*ps++) + zero;
 
850
                        dreg = (double) *pf++;
 
851
                        if (dreg<mnvd)
 
852
                           mnvd = dreg;
 
853
                        else if (mxvd<dreg)
 
854
                           mxvd = dreg;
 
855
                        }
 
856
                     }
 
857
                  else
 
858
                     {
 
859
                     for (k=0; k<dn; k++)
 
860
                        {
 
861
                        *pf = (*ps++) + zero;
 
862
                        dreg = (double) *pf++;
 
863
                        if (dreg<mnvd)
 
864
                           mnvd = dreg;
 
865
                        else if (mxvd<dreg)
 
866
                           mxvd = dreg;
 
867
                        }
 
868
                     }
 
869
                  }
 
870
               }
 
871
            else
 
872
               {
 
873
               if (bf)
 
874
                  {
 
875
                  if (ifac == 0)
 
876
                     {
 
877
                     for (k=0; k<dn; k++, pf++, ps++) 
 
878
                        {
 
879
                        if (*ps==s) 
 
880
                           toNULLFLOAT(*pf);
 
881
                        else
 
882
                           *pf = fac*(*ps) + zero;
 
883
                        }
 
884
                     }
 
885
                  else
 
886
                     {
 
887
                     for (k=0; k<dn; k++, pf++, ps++) 
 
888
                        {
 
889
                        if (*ps==s) 
 
890
                           toNULLFLOAT(*pf);
 
891
                        else
 
892
                           *pf = *ps + zero;
 
893
                        }
 
894
                     }
 
895
                  }
 
896
               else
 
897
                  {
 
898
                  if (ifac == 0)
 
899
                     {
 
900
                     for (k=0; k<dn; k++)
 
901
                           *pf++ = fac*(*ps++) + zero;
 
902
                     }
 
903
                  else
 
904
                     {
 
905
                     for (k=0; k<dn; k++) *pf++ = (*ps++) + zero;
 
906
                     }
 
907
                  }
 
908
               }
 
909
 
 
910
            if (Midas_flag != 0)
 
911
               MyPut(-32,felem,dn,(char *) buf.f);
 
912
            else
 
913
               SCFPUT(mfd,felem,dn,(char *) buf.f);
 
914
            }
 
915
 
 
916
         /* option = convert data to float (unsingend int and int are equal) */
 
917
         /* format for minmax is float */
 
918
 
 
919
         else if (fmt=='F') 
 
920
            {   
 
921
            pf = buf.f; s = bfdef->blank, ps = p.s;
 
922
            if (fmm)
 
923
               {                        /* dsf = -32 */
 
924
               if (felem==1) mnvd = mxvd = (float) *ps;
 
925
               if (bf)
 
926
                  {
 
927
                  for (k=0; k<dn; k++, pf++, ps++)
 
928
                     {
 
929
                     if (*ps == s)
 
930
                        toNULLFLOAT(*pf);
 
931
                     else
 
932
                        {
 
933
                        *pf = (float) *ps;
 
934
                        dreg = (double) *pf;
 
935
                        if (dreg<mnvd)
 
936
                           mnvd = dreg;
 
937
                        else if (mxvd<dreg)
 
938
                           mxvd = dreg;
 
939
                        }
 
940
                     }
 
941
                  }
 
942
               else
 
943
                  {
 
944
                  for (k=0; k<dn; k++) 
 
945
                     {
 
946
                     *pf = (float) *ps++;
 
947
                     dreg = (double) *pf++;
 
948
                     if (dreg<mnvd)
 
949
                        mnvd = dreg;
 
950
                     else if (mxvd<dreg)
 
951
                        mxvd = dreg;
 
952
                     }
 
953
                  }
 
954
               }
 
955
            else
 
956
               {
 
957
               if (bf)
 
958
                  {
 
959
                  for (k=0; k<dn; k++, pf++, ps++)
 
960
                     {
 
961
                     if (*ps==s)
 
962
                        toNULLFLOAT(*pf);
 
963
                     else
 
964
                        *pf = (float) *ps;
 
965
                     }
 
966
                  }
 
967
               else
 
968
                  {
 
969
                  for (k=0; k<dn; k++) *pf++ = (float) *ps++;
 
970
                  }
 
971
               }
 
972
 
 
973
            if (Midas_flag != 0)
 
974
               MyPut(-32,felem,dn,(char *) buf.f);
 
975
            else
 
976
               SCFPUT(mfd,felem,dn,(char *) buf.f);
 
977
            }
 
978
 
 
979
         /* option = take data as is (unsigned int and int are different) */
 
980
         /* format for minmax is same as original format */
 
981
 
 
982
         else 
 
983
            {           
 
984
            if (dfmt == (-16))
 
985
               {
 
986
               pu = p.u; u = bfdef->blank; ps = p.s;
 
987
               if (fmm)
 
988
                  {                             /* here dsf = dfmt */
 
989
                  if (felem==1) mnvi = mxvi = *ps + 32768;
 
990
                  if (bf)
 
991
                     {
 
992
                     for (k=0; k<dn; k++, pu++)
 
993
                        {
 
994
                        ireg = (int) *ps++;
 
995
                        if (ireg == u)
 
996
                           toNULLSHORT(*pu);
 
997
                        else
 
998
                           {
 
999
                           ireg += 32768;
 
1000
                           if (ireg < mnvi)
 
1001
                              mnvi = ireg;
 
1002
                           else if (mxvi<ireg)
 
1003
                              mxvi = ireg;
 
1004
                           *pu = ireg;
 
1005
                           }
 
1006
                        }
 
1007
                     }
 
1008
                  else
 
1009
                     {
 
1010
                     for (k=0; k<dn; k++) 
 
1011
                        {
 
1012
                        ireg = (int) *ps++;
 
1013
                        ireg += 32768;
 
1014
                        if (ireg < mnvi)
 
1015
                           mnvi = ireg;
 
1016
                        else if (mxvi<ireg)
 
1017
                           mxvi = ireg;
 
1018
                        *pu++ = ireg;
 
1019
                        }
 
1020
                     }
 
1021
                  }
 
1022
               else
 
1023
                  {
 
1024
                  if (bf)
 
1025
                     {
 
1026
                     for (k=0; k<dn; k++, pu++, ps++) 
 
1027
                        {
 
1028
                        if (*ps == u) 
 
1029
                           toNULLSHORT(*pu);
 
1030
                        else
 
1031
                           *pu = *ps + 32768;
 
1032
                        }
 
1033
                     }
 
1034
                  else
 
1035
                     {
 
1036
                     for (k=0; k<dn; k++) *pu++ = *ps++ + 32768;
 
1037
                     }
 
1038
                  }
 
1039
               if (Midas_flag != 0)
 
1040
                  MyPut(dfmt,felem,dn,(char *) p.u);
 
1041
               else
 
1042
                  SCFPUT(mfd,felem,dn,(char *) p.u);
 
1043
               }
 
1044
            else                /* dfmt = 16 */
 
1045
               {
 
1046
               ps = p.s; s = bfdef->blank;
 
1047
               if (fmm)
 
1048
                  {                          /* here dsf = dfmt */
 
1049
                  if (felem==1) mnvi = mxvi = *ps;
 
1050
                  if (bf)
 
1051
                     {
 
1052
                     for (k=0; k<dn; k++, ps++)
 
1053
                        {
 
1054
                        ireg = (int) *ps;
 
1055
                        if (ireg == s) 
 
1056
                           toNULLSHORT(*ps);
 
1057
                        else
 
1058
                           {
 
1059
                           if (ireg<mnvi)
 
1060
                              mnvi = ireg;
 
1061
                           else if (mxvi<ireg)
 
1062
                              mxvi = ireg;
 
1063
                           }
 
1064
                        }
 
1065
                     }
 
1066
                  else
 
1067
                     {
 
1068
                     for (k=0; k<dn; k++)
 
1069
                        {
 
1070
                        ireg = (int) *ps++;
 
1071
                        if (ireg<mnvi)
 
1072
                           mnvi = ireg;
 
1073
                        else if (mxvi<ireg)
 
1074
                           mxvi = ireg;
 
1075
                        }
 
1076
                     }
 
1077
                  }
 
1078
               else
 
1079
                  {
 
1080
                  if (bf) 
 
1081
                     {
 
1082
                     for (k=0; k<dn; k++, ps++)
 
1083
                        {
 
1084
                        if (*ps == s) toNULLSHORT(*ps);
 
1085
                        }
 
1086
                     }
 
1087
                  }
 
1088
 
 
1089
               if (Midas_flag != 0)
 
1090
                  MyPut(dfmt,felem,dn,(char *) p.s);
 
1091
               else
 
1092
                  SCFPUT(mfd,felem,dn,(char *) p.s);
 
1093
               }
 
1094
            }
 
1095
 
 
1096
         felem += dn;
 
1097
         if (!dno) 
 
1098
            {
 
1099
            gno++;
 
1100
            pno = bfdef->pcount; pcnt = 0;
 
1101
            pp = bfdef->parm;
 
1102
            dno = ndata;
 
1103
            }
 
1104
         }
 
1105
 
 
1106
      } while (npix && gno<bfdef->gcount);
 
1107
   }
 
1108
 
 
1109
 
 
1110
if (fmm) 
 
1111
   {                         /* save computed min/max values  */
 
1112
   if (dsf > -32) 
 
1113
      {                         /* handle also UI2 here */
 
1114
      cuts[0] = (float)mnvi; cuts[1] = (float)mxvi;
 
1115
      }
 
1116
   else 
 
1117
      {
 
1118
      if (MAXFLOAT<mnvd) 
 
1119
         mnvd = MAXFLOAT;
 
1120
      else if (mnvd<MINFLOAT) 
 
1121
         mnvd = MINFLOAT;
 
1122
      if (MAXFLOAT<mxvd) 
 
1123
         mxvd = MAXFLOAT;
 
1124
      else if (mxvd<MINFLOAT) 
 
1125
         mxvd = MINFLOAT;
 
1126
      cuts[0] = (float)mnvd; cuts[1] = (float)mxvd;
 
1127
      }
 
1128
    SCDWRR(mfd,"LHCUTS",cuts,3,2,unit);
 
1129
 
 
1130
   }
 
1131
 
 
1132
 
 
1133
if (Midas_flag == 0)
 
1134
   {
 
1135
   if (0<=mfd) SCFCLO(mfd);                       /* close data files  */
 
1136
   if (0<=mfdt) { TCSINI(mfdt); TCTCLO(mfdt); }
 
1137
   mfd = -1; mfdt = -1;
 
1138
   }
 
1139
 
 
1140
return 0;
 
1141
}
 
1142
 
 
1143