← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/plot/src/plotvec.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/plot/src/plotvec.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1993-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
.IDENTifer   PLOTVEC
 
30
.AUTHOR      R.M. van Hees IPG-ESO Garching
 
31
.KEYWORDS    Graphics, bulk data frame
 
32
.LANGUAGE    C
 
33
.PURPOSE     Plots vectors (polarisation map) from 2 two-dimensional FRAME
 
34
  in/output: IN_A/C/1/60   = input frame with intensisty information
 
35
             IN_B/C/1/60   = input frame with position angle information
 
36
             INPUTC/C/1/72 = coord_string (default: whole frame)
 
37
             INPUTR/R/1/3  = scale of vector [units/mm]
 
38
                             position angle range (default: 0,360)
 
39
             INPUTI/I/1/2  = smooting parameter (default: 1)
 
40
                             head (default: 1, the vector get an arrow)
 
41
 
 
42
.COMMENTS    none
 
43
.ENVIRONment MIDAS
 
44
             #include <midas_def.h>     Prototypes for MIDAS interfaces
 
45
             #include <plot_def.h>      General symbols for Plot routines
 
46
 
 
47
.VERSION     1.1     25-Nov-1993   Removed restrictions on frame size, RvH
 
48
             1.0     10-Sep-1993   FORTRAN --> C        RvH
 
49
 
 
50
 090710         last modif
 
51
------------------------------------------------------------*/
 
52
/*
 
53
 * Define _POSIX_SOURCE to indicate
 
54
 * that this is a POSIX program
 
55
 */
 
56
#define _POSIX_SOURCE 1
 
57
 
 
58
/*
 
59
 * definition of the used functions
 
60
 */
 
61
#include <stdio.h>
 
62
#include <string.h>
 
63
#include <stdlib.h>
 
64
#include <math.h>
 
65
 
 
66
/*
 
67
 * define some macros and constants
 
68
 */
 
69
#include <midas_def.h>
 
70
#include <plot_def.h>
 
71
 
 
72
#define MAXLEV          50
 
73
#define MAXPIX          512    /* max frame dimension (X,Y) accessed at once */
 
74
#define MAXSIZ          (MAXPIX * MAXPIX)
 
75
 
 
76
/*
 
77
 * here starts the code of the function
 
78
 */
 
79
int main()
 
80
{
 
81
register int ii;
 
82
int     actvals, chunks, head, imfA, imfB, knul, last_row, naxis, 
 
83
        size, stat, unit, ndum[PLDIM2], npixA[PLDIM2], npixB[PLDIM2], 
 
84
        sublo[PLDIM2], subhi[PLDIM2];
 
85
 
 
86
float   amin, amax, *p_imgA, *p_imgB, scar, nrpix[2], ranp[4], area[4], 
 
87
        image[4], wcfram[12];
 
88
 
 
89
double  startA[PLDIM2], stepA[PLDIM2], startB[PLDIM2], stepB[PLDIM2];
 
90
 
 
91
char    cmnd[21], nameA[61], nameB[61], cunitA[49], identA[33], identB[33], 
 
92
        input[73], *label[4];
 
93
/*
 
94
 * initialised variables
 
95
 */
 
96
char  *err_start = "*** FATAL: START descriptor of the frames differ",
 
97
      *err_step  = "*** FATAL: STEP descriptor of the frames differ",
 
98
      *err_npix  = "*** FATAL: NPIX descriptor of the frames differ",
 
99
      *err_1dim  = "*** FATAL: the frames have only one dimension",
 
100
      *err_coord = "*** FATAL: invalid coordinate input ...",
 
101
      *err_xcuts = "*** FATAL: range in x has no overlap with current graph abscissa - NO PLOT",
 
102
      *err_ycuts = "*** FATAL: range in y has no overlap with current graph abscissa - NO PLOT";
 
103
 
 
104
static char *axis[PLDIM2] = { "MANU", "MANU" };
 
105
 
 
106
int     access =  0,                      /* parameter for PCOPEN: plot mode */
 
107
        plmode = -1,                   /* plot mode taken from keyword PMODE */
 
108
        ismoot =  1;                      /* Default no smooting of the data */
 
109
/*
 
110
 * allocate memory for different character pointers and initialise a few
 
111
 */
 
112
for ( ii = 0; ii < 4; ii++ ) label[ii] = osmmget(81);
 
113
 
 
114
(void) strcpy( label[0], "Position (" );
 
115
(void) strcpy( label[1], "Position (" );
 
116
 
 
117
/*
 
118
 * start of executable code
 
119
 */
 
120
(void) SCSPRO( "PLTVEC" );                   /*contact with the MIDAS monitor*/
 
121
(void) SCKGETC( "MID$CMND", 1, 20, &actvals, cmnd );
 
122
if ( *cmnd == 'O' ) access = 1; 
 
123
 
 
124
/*
 
125
 * find file name and size information of first frame
 
126
 */
 
127
(void) SCKGETC( "IN_A", 1, 60, &actvals, nameA );
 
128
(void) SCFOPN( nameA, D_R4_FORMAT, 0, F_IMA_TYPE, &imfA );
 
129
(void) SCDRDI( imfA, "NAXIS", 1, 1, &actvals, &naxis, &unit, &knul );
 
130
(void) SCDRDI( imfA, "NPIX" , 1, PLDIM2, &actvals, npixA , &unit, &knul );
 
131
if ( naxis < 2 || (npixA[0] == 1 || npixA[1] == 1) ) SCETER( 1, err_1dim );
 
132
 
 
133
/*
 
134
 * find file name and size information of second frame
 
135
 */
 
136
(void) SCKGETC( "IN_B", 1, 60, &actvals, nameB );
 
137
(void) SCFOPN( nameB, D_R4_FORMAT, 0, F_IMA_TYPE, &imfB );
 
138
(void) SCDRDI( imfB, "NPIX" , 1, PLDIM2, &actvals, npixB , &unit, &knul );
 
139
if ( npixA[0] != npixB[0] || npixA[1] != npixB[1] ) SCETER( 2, err_npix );
 
140
 
 
141
/*
 
142
 * read descriptor of first frame
 
143
 */
 
144
(void) SCDRDD( imfA, "START", 1, PLDIM2, &actvals, startA, &unit, &knul );
 
145
(void) SCDRDD( imfA, "STEP" , 1, PLDIM2, &actvals, stepA , &unit, &knul );
 
146
(void) SCDGETC( imfA, "IDENT", 1, 32, &actvals, identA );
 
147
(void) SCDGETC( imfA, "CUNIT", 1, 48, &actvals, cunitA );
 
148
 
 
149
/*
 
150
 * read descriptor of second frame
 
151
 */
 
152
(void) SCDRDD( imfB, "START", 1, PLDIM2, &actvals, startB, &unit, &knul );
 
153
(void) SCDRDD( imfB, "STEP" , 1, PLDIM2, &actvals, stepB , &unit, &knul );
 
154
(void) SCDGETC( imfB, "IDENT", 1, 32, &actvals, identB );
 
155
 
 
156
/*
 
157
 * do both frames overlap?
 
158
 */
 
159
if ( startA[0] != startB[0] || startA[1] != startB[1] ) SCETER( 3, err_start );
 
160
if ( stepA[0] != stepB[0] || stepA[1] != stepB[1] ) SCETER( 4, err_step );
 
161
 
 
162
/*
 
163
 * Get the manual setting for the axes
 
164
 */
 
165
PCKRDR( "XAXIS", 4, &actvals, wcfram );
 
166
PCKRDR( "YAXIS", 4, &actvals, wcfram+FOR_Y );
 
167
 
 
168
/*
 
169
 * read window coordinates and take action
 
170
 */
 
171
(void) SCKGETC( "INPUTC", 1, 60, &actvals, input );
 
172
 
 
173
if ( *input == 'm' || *input == 'M' )                      /* manual scaling */
 
174
   { BOXWTP( wcfram,       npixA[0], startA[0], stepA[0], image );
 
175
     BOXWTP( wcfram+FOR_Y, npixA[1], startA[1], stepA[1], image+2 );
 
176
   }
 
177
else
 
178
   { if ( *input == 'c' || *input == 'C' )                  /* display input */
 
179
        { (void) SCKRDR( "OUTPUTR", 10, 1, &actvals, image, &unit, &knul);
 
180
          (void) SCKRDR( "OUTPUTR", 11, 1, &actvals, image+2, &unit, &knul); 
 
181
          (void) SCKRDR( "OUTPUTR", 15, 1, &actvals, image+1, &unit, &knul);
 
182
          (void) SCKRDR( "OUTPUTR", 16, 1, &actvals, image+3, &unit, &knul);
 
183
        }
 
184
     else                                               /* automatic scaling */
 
185
        { stat = Convcoo(1,imfA,input,PLDIM2,ndum,sublo,subhi);
 
186
          if ( stat != ERR_NORMAL ) SCETER( 5, err_coord );
 
187
          image[0] = sublo[0] + 1;
 
188
          image[1] = subhi[0] + 1;
 
189
          image[2] = sublo[1] + 1;
 
190
          image[3] = subhi[1] + 1;
 
191
        }
 
192
   }
 
193
BOXPTW( image,   npixA[0], startA[0], stepA[0], area );
 
194
BOXPTW( image+2, npixA[1], startA[1], stepA[1], area+2 );
 
195
 
 
196
PCKWRR( "PIXEL", 4, image );
 
197
 
 
198
if ( access == 0 )
 
199
   { 
 
200
/*
 
201
 * get size of frame along X-axis
 
202
 */
 
203
     if ( fabs( *wcfram ) < PLT_EPS && fabs( *(wcfram+1) ) < PLT_EPS )
 
204
        { wcfram[0] = area[2];
 
205
          wcfram[1] = area[3]; 
 
206
          wcfram[2] = wcfram[3] = 0.0;
 
207
        }
 
208
/*
 
209
 * get size of frame along Y-axis
 
210
 */
 
211
     if ( fabs( *(wcfram+FOR_Y) ) < PLT_EPS 
 
212
          && fabs( *(wcfram+FOR_Y+1) ) < PLT_EPS )
 
213
        { axis[1] = "AUTO";
 
214
          wcfram[FOR_Y]   = area[2];
 
215
          wcfram[FOR_Y+1] = area[3]; 
 
216
          wcfram[FOR_Y+2] = wcfram[FOR_Y+3] = 0.0;
 
217
        }
 
218
     GETFRM( axis[0], wcfram );
 
219
     GETFRM( axis[1], wcfram + FOR_Y );
 
220
     PCKWRR( "XWNDL", 4, wcfram );
 
221
     PCKWRR( "YWNDL", 4, wcfram+FOR_Y );
 
222
   }
 
223
else                                                         /* overplot mode*/
 
224
   { PCKRDR( "XWNDL", 4, &actvals, wcfram );
 
225
     PCKRDR( "YWNDL", 4, &actvals, wcfram+FOR_Y );                       
 
226
/*
 
227
 * does overplot data  fall within plotted frame? 
 
228
 */
 
229
     amin = MYMIN( *wcfram, *(wcfram + 1) );
 
230
     amax = MYMAX( *wcfram, *(wcfram + 1) );
 
231
     if ( ( MYMAX( area[0], area[1] ) < amin ) ||
 
232
          ( MYMIN( area[0], area[1] ) > amax ) )
 
233
        SCETER( 6, err_xcuts );
 
234
 
 
235
     amin = MYMIN( *(wcfram + FOR_Y), *(wcfram + FOR_Y + 1) );
 
236
     amax = MYMAX( *(wcfram + FOR_Y), *(wcfram + FOR_Y + 1) );
 
237
     if ( ( MYMAX( area[2], area[3] ) < amin ) ||
 
238
          ( MYMIN( area[2], area[3] ) > amax ) )
 
239
        SCETER( 7, err_ycuts );
 
240
   }
 
241
/*
 
242
 * get intensity
 
243
 */
 
244
(void) SCDRDR( imfA, "LHCUTS", 1, 4, &actvals, wcfram+FOR_Z, &unit, &knul);
 
245
(void) SCKRDR( "INPUTR", 2, 4, &actvals, ranp, &unit, &knul);
 
246
if ( fabs(ranp[0]) > PLT_EPS || fabs(ranp[1]) > PLT_EPS )
 
247
   { wcfram[FOR_Z]  = ranp[0];
 
248
     wcfram[FOR_Z+1]= ranp[1];
 
249
   }
 
250
else
 
251
   {
 
252
   if ( wcfram[FOR_Z] == wcfram[FOR_Z+1] )
 
253
      { wcfram[FOR_Z] = wcfram[FOR_Z+2];
 
254
        wcfram[FOR_Z+1] = wcfram[FOR_Z+3];
 
255
      }
 
256
    }
 
257
ranp[0] = wcfram[FOR_Z];
 
258
ranp[1] = wcfram[FOR_Z+1];
 
259
PCKWRR( "ZWNDL", 2, wcfram+FOR_Z);
 
260
(void) SCKRDR( "INPUTR", 1, 1, &actvals, &scar, &unit, &knul);
 
261
/*
 
262
 * the intensity scale
 
263
 */
 
264
if ( scar == 0.0 )
 
265
   { if ( wcfram[FOR_Z] == wcfram[FOR_Z+1] )
 
266
        scar = wcfram[FOR_Z] / 10.0;
 
267
     else
 
268
        scar = ( wcfram[FOR_Z+1] - wcfram[FOR_Z] ) / 10.0;
 
269
   }
 
270
/*
 
271
 * get the other input parameters
 
272
 */
 
273
(void) SCKRDI( "INPUTI", 1, 1, &actvals, &ismoot, &unit, &knul);
 
274
(void) SCKRDI( "INPUTI", 2, 1, &actvals, &head, &unit, &knul);
 
275
 
 
276
/*
 
277
 * setup graphic device according to MIDAS settings
 
278
 */
 
279
PCOPEN( " ", " ", access, &plmode );
 
280
 
 
281
/*
 
282
 * determine number of pixels and number of chunks
 
283
 */
 
284
nrpix[0] = (int) fabs( image[1] - image[0] ) + 1;
 
285
nrpix[1] = (int) fabs( image[3] - image[2] ) + 1;
 
286
last_row = MYMAX( image[2], image[3]);
 
287
chunks   = (int) ceil( (double) nrpix[0] * nrpix[1] / MAXSIZ );
 
288
nrpix[1] = (int) ceil( (double) nrpix[1] / chunks );
 
289
 
 
290
/*
 
291
 * allocate virtual memory and scratch space
 
292
 */
 
293
size  = nrpix[0] * nrpix[1];
 
294
p_imgA = (float *) osmmget( size * sizeof( float ));
 
295
p_imgB = (float *) osmmget( size * sizeof( float ));
 
296
 
 
297
for ( ii = 0; ii < chunks; ii++ )
 
298
    { if ( image[3] > image[2] )                        /* size of chunk p.c.*/
 
299
         { if ( ii > 0 ) image[2] += nrpix[1] - 1.0;
 
300
           image[3] = MYMIN( last_row, image[2] + nrpix[1] - 1 );
 
301
         }
 
302
      else
 
303
         { if ( ii > 0 ) image[3] += nrpix[1] - 1.0;
 
304
           image[2] = MYMIN( last_row, image[3] + nrpix[1] - 1 );
 
305
         }
 
306
/*
 
307
 * get size of chunk in w.c.
 
308
 */
 
309
      BOXPTW( image+2, npixA[1], startA[1], stepA[1], area+2 );      
 
310
/*
 
311
 * extract chunck from original frame
 
312
 */
 
313
      GETDAT( imfA, MAXSIZ, npixA, image, ismoot, p_imgA );
 
314
      GETDAT( imfB, MAXSIZ, npixB, image, ismoot, p_imgB );
 
315
/*
 
316
 * do it!!
 
317
 */
 
318
      PLVEC( p_imgA, p_imgB, image, area, stepA, scar, ranp, head );
 
319
/*
 
320
 * updating for the next loop
 
321
 */
 
322
      nrpix[1] = MYMIN( last_row - ii * nrpix[1], nrpix[1] );
 
323
      size  = nrpix[0] * nrpix[1];
 
324
    }
 
325
(void) SCFCLO( imfA );
 
326
(void) SCFCLO( imfB );
 
327
 
 
328
/*
 
329
 * draw the axes and the label
 
330
 */
 
331
if ( plmode >= 0 && access == 0 )
 
332
   { if ( strlen( cunitA ) > (size_t) 32 ) 
 
333
        { (void) strncat( label[0], cunitA+32, 16 );
 
334
          *(cunitA+32) = '\0';
 
335
        }
 
336
     if ( strlen( cunitA ) > (size_t) 16 ) 
 
337
        (void) strcat( label[1], cunitA+16 );
 
338
 
 
339
     for ( ii = 0; ii < PLDIM2; ii++ )
 
340
         { (void) strcat( label[ii], ")" );
 
341
           LABSTR( label[ii] );
 
342
         }
 
343
/*
 
344
 * plot axes and labels
 
345
 */
 
346
     PCFRAM( wcfram, wcfram+FOR_Y, label[0], label[1] );
 
347
 
 
348
     if ( plmode == 1 )
 
349
        { (void) sprintf( label[2], "Frames: %s + %s", nameA, nameB );
 
350
          (void) sprintf( label[3], "Ident: %s + %s", identA, identB );
 
351
          PLIDEN( plmode, label[2], label[3] );
 
352
        }
 
353
     else if ( plmode == 2 )
 
354
        PLVECI( plmode, nameA, identA, nameB, identB, scar, ranp );
 
355
   }
 
356
/*
 
357
 * close plot file and terminate graphic operations
 
358
 */
 
359
PCCLOS();
 
360
 
 
361
return SCSEPI();
 
362
}
 
363