1 /* Copyright (C) 1988-2005 by Brian Doty and the Institute
2 of Global Environment and Society (IGES).
3
4 See file COPYRIGHT for more information. */
5
6 /* Authored by B. Doty */
7
8 #ifdef HAVE_CONFIG_H
9 #include <config.h>
10
11 /* If autoconfed, only include malloc.h when it's presen */
12 #ifdef HAVE_MALLOC_H
13 #include <malloc.h>
14 #endif
15
16 #else /* undef HAVE_CONFIG_H */
17
18 #include <malloc.h>
19
20 #endif /* HAVE_CONFIG_H */
21
22 #include <stdio.h>
23 #include <math.h>
24 #include <string.h>
25 #include <stdlib.h>
26 #include <time.h>
27 #include <errno.h>
28 #include "grads.h"
29 #include "gx.h"
30 #include "wx.h"
31
32 void gatmlb (struct gacmn *); /*mf 970727 -- time label ---mf*/
33
34 static char pout[256]; /* Build error msgs here */
35 static struct mapprj mpj; /* Common map projection structure */
36
37 static float wxymin,wxymax; /* wx symbol limits */
38
39 /* Current I and J axis translation routines for grid to absolute
40 for use by gaconv */
41
42 /* A note: the conversion routines normally handle grid
43 coordinates relative to the file (for convenience).
44 But the graphics routines normally use grid coordinates
45 with respect to the grid being worked with (ie, values
46 1 to n). Thus the gaconv routine handles this translation
47 with the ioffset and joffset values. */
48
49 static float (*iconv) (float *, float);
50 static float (*jconv) (float *, float);
51 static float *ivars, *jvars;
52 static float ioffset, joffset;
53 static float idiv, jdiv;
54
gagx(struct gacmn * pcm)55 void gagx (struct gacmn *pcm) {
56 gxstrt (pcm->xsiz,pcm->ysiz,pcm->batflg,pcm->hbufsz);
57 pcm->pass = 0;
58 pcm->ccolor = -9;
59 pcm->cint = 0.0;
60 pcm->cstyle = -9;
61 pcm->cthick = 3;
62 pcm->shdcnt = 0;
63 pcm->lastgx = 0;
64 pcm->xdim = -1;
65 pcm->ydim = -1;
66 pcm->xgr2ab = NULL;
67 pcm->ygr2ab = NULL;
68 pcm->xab2gr = NULL;
69 pcm->yab2gr = NULL;
70 }
71
72 int rcols[13] = {9,14,4,11,5,13,3,10,7,12,8,2,6};
73
74 /* Figure out which graphics routine to call. Use the first
75 grid hung off gacmn to determine whether we are doing
76 a 0-D, 1-D, or 2-D output. */
77
gaplot(struct gacmn * pcm)78 void gaplot (struct gacmn *pcm) {
79 struct gagrid *pgr;
80 struct gastn *stn;
81 int proj;
82 int rc;
83
84 pcm->relnum = pcm->numgrd;
85 proj = pcm->mproj;
86 if (pcm->mproj>1 && (pcm->xflip || pcm->yflip)) pcm->mproj = 1;
87
88 /* if gxout stat, then doesn't matter what type grid */
89
90 if (pcm->gout0==1) {
91 gastts(pcm);
92 } else if (pcm->gout0==2) {
93 gadprnt(pcm);
94 } else if (pcm->gout0==3) {
95 gaoutgds(pcm);
96 } else {
97
98 /* If statflg, produce gxout-stat output for all displays */
99
100 if (pcm->statflg) gastts(pcm);
101
102 /* output is a grid */
103
104 if (pcm->type[0] == 1) {
105 pgr = pcm->result[0].pgr;
106 if (pgr->idim==-1) { /* 0-D */
107 if (pcm->gout2a == 7 ) {
108 gafwrt (pcm);
109 #if USELATS == 1
110 } else if (pcm->gout2a == 20 ) {
111 rc=galats(pcm,3,0);
112 } else if (pcm->gout2a == 21 ) {
113 rc=galats(pcm,5,0);
114 #endif
115 } else {
116 sprintf (pout,"Result value = %g \n",pgr->rmin);
117 gaprnt (2,pout);
118 }
119 }
120 else if (pgr->jdim==-1) { /* 1-D */
121 if (pcm->gout2a==7) {
122 gafwrt (pcm);
123 #if USELATS == 1
124 } else if (pcm->gout2a==20) {
125 rc=galats(pcm,3,0);
126 } else if (pcm->gout2a==21) {
127 rc=galats(pcm,5,0);
128 #endif
129 } else if (pcm->gout2b==5 && pcm->numgrd>1) gascat(pcm);
130 else if (pcm->gout1==1) gagrph(pcm,0);
131 else if (pcm->gout1==2) gagrph(pcm,1);
132 else if (pcm->gout1==3) gagrph(pcm,2);
133 else galfil(pcm);
134 }
135 else { /* 2-D */
136 if (pcm->numgrd==1) {
137 if (pcm->gout2a == 1 ) {
138 gacntr (pcm,0);
139 } else if (pcm->gout2a == 2 ) {
140 gacntr (pcm,1);
141 } else if (pcm->gout2a == 3 ) {
142 gaplvl (pcm);
143 } else if (pcm->gout2a == 6 ) {
144 gafgrd (pcm);
145 } else if (pcm->gout2a == 7 ) {
146 gafwrt (pcm);
147 #if USELATS == 1
148 } else if (pcm->gout2a == 20 ) {
149 rc=galats(pcm,3,0);
150 } else if (pcm->gout2a == 21 ) {
151 rc=galats(pcm,5,0);
152 #endif
153 } else {
154 gacntr (pcm,2);
155 }
156 } else {
157 if (pcm->gout2b == 3 ) {
158 gaplvl (pcm);
159 } else if (pcm->gout2b == 4 ) {
160 gavect (pcm,0);
161 } else if (pcm->gout2b == 5 ) {
162 gascat (pcm);
163 } else if (pcm->gout2b == 8 ) {
164 gastrm (pcm);
165 } else {
166 gavect (pcm,1);
167 }
168 }
169 }
170
171 /* Output is station data */
172
173 } else {
174 stn = pcm->result[0].stn;
175 if (stn->idim==0 && stn->jdim==1) {
176 if (pcm->goutstn==1 || pcm->goutstn==2 || pcm->goutstn==6)
177 gapstn (pcm);
178 else if (pcm->goutstn==3) gafstn(pcm);
179 else if (pcm->goutstn==4) gapmdl (pcm);
180 else if (pcm->goutstn==7) gasmrk (pcm);
181 else gawsym (pcm);
182 }
183 else if (stn->idim==2 && stn->jdim == -1) gapprf (pcm);
184 else if (stn->idim==3 && stn->jdim == -1) gatser (pcm);
185 else {
186 gaprnt (0,"Invalid station data dimension environment\n");
187 }
188 }
189 }
190 pcm->mproj = proj;
191 }
192
gawgdsval(FILE * outfile,float * val)193 void gawgdsval(FILE* outfile, float *val) {
194 if (BYTEORDER != 1) { /* always write big endian for the GDS */
195 gabswp(val, 1);
196 }
197 fwrite(val, sizeof(float), 1, outfile);
198 }
199
gawgdstime(FILE * outfile,double * val)200 void gawgdstime(FILE* outfile, double *val) {
201 float tmp;
202 sprintf(pout, "pre-byteswapped time: %f", *val); gaprnt(0, pout);
203 if (BYTEORDER != 1) { /* always write big endian for the GDS */
204 ganbswp(val, sizeof(double));
205 }
206 sprintf(pout, "byteswapped time: %f", *val); gaprnt(0, pout);
207 fwrite(val, sizeof(double), 1, outfile);
208 }
209
210
211 /* Writes station data out to a file as a DODS sequence */
gaoutgds(struct gacmn * pcm)212 void gaoutgds (struct gacmn *pcm) {
213
214 const char startrec[4] = {0x5a, 0x00, 0x00, 0x00};
215 const char stnidlen[4] = {0x00, 0x00, 0x00, 0x08};
216 const char endrec[4] = {0xa5, 0x00, 0x00, 0x00};
217
218 int rc, i, numvars, retval, levelstart, numreps;
219 char *sendstnid, *sendlat, *sendlon, *sendlev, *sendtime, *senddep, *sendind;
220 struct garpt **currpt;
221 struct garpt *ref, *levelref;
222 int *varlevels;
223 double coardstime;
224 FILE *outfile;
225 char *sendoptions;
226 float outFloat;
227
228 if (pcm->wgds->fname == NULL) {
229 gaprnt(0, "error: no file specified (use \"set writegds\").\n");
230 return;
231 }
232 outfile = fopen(pcm->wgds->fname, "ab");
233 if (outfile == NULL) {
234 gaprnt(0, "error: WRITEGDS unable to open ");
235 gaprnt(0, pcm->wgds->fname);
236 gaprnt(0, " for write\n");
237 return;
238 }
239 gaprnt(0, "got options and opened file\n");
240
241 if (pcm->wgds->opts == NULL) {
242 gaprnt(0, "No options specified. Defaulting to full output (\"sxyztdi\").\n");
243 sendoptions = "sxyztdi";
244 } else {
245 sendoptions = pcm->wgds->opts;
246 }
247
248 if (strchr(sendoptions, 'f')) {
249 /* Write a DODS End-Of-Sequence marker to indicate to the client
250 * that there is no more data. This is separate from gaoutgds() so
251 * that time loops can be written as a single sequence.
252 */
253 gaprnt(0, "Finishing sequence:\n");
254 gaprnt(0, "EOS\n");
255 fwrite(endrec, sizeof(char), 4, outfile);
256 fclose(outfile);
257 return;
258 }
259
260 sendstnid = strchr(sendoptions, 's');
261 sendlon = strchr(sendoptions, 'x');
262 sendlat = strchr(sendoptions, 'y');
263 sendlev = strchr(sendoptions, 'z');
264 sendtime = strchr(sendoptions, 't');
265 senddep = strchr(sendoptions, 'd');
266 sendind = strchr(sendoptions, 'i');
267
268 numvars = pcm->numgrd;
269 currpt = (struct garpt **)malloc(sizeof(struct garpt *) * numvars);
270 varlevels = (int *)malloc(sizeof(int) * numvars);
271 if (currpt == NULL || varlevels == NULL) {
272 gaprnt(0, "error: memory allocation failed\n");
273 fclose(outfile);
274 return;
275 }
276 for (i = 0; i < numvars; i++) {
277 varlevels[i] = pcm->result[i].stn->pvar->levels;
278 currpt[i] = pcm->result[i].stn->rpt;
279 }
280
281 /* set levelstart to index of first level-dependent variable.
282 * if none is found, levelstart = numvars (this is used below) */
283 levelstart = 0;
284 while (varlevels[levelstart] == 0 && levelstart < numvars) {
285 levelstart++;
286 }
287
288 retval = 0;
289 numreps = 0;
290
291 /* write reports */
292 while (currpt[0]) {
293 ref = currpt[0];
294 coardstime = ref->tim; /* change to meaningful conversion */
295
296 gaprnt(0, "SOI ");
297 fwrite(startrec, sizeof(char), 4, outfile);
298
299 gaprnt(0, ">>\t");
300
301 if (sendstnid) {
302 sprintf(pout, "stnid: %.8s ", ref->stid); gaprnt(0, pout);
303 fwrite(stnidlen, sizeof(char), 4, outfile);
304 fwrite(&(ref->stid), sizeof(char), 8, outfile);
305 }
306 if (sendlon) {
307 sprintf(pout, "lon: %f ", ref->lon); gaprnt(0, pout);
308 outFloat = ref->lon;
309 gawgdsval(outfile, &outFloat);
310 }
311 if (sendlat) {
312 sprintf(pout, "lat: %f ", ref->lat); gaprnt(0, pout);
313 outFloat = ref->lat;
314 gawgdsval(outfile, &outFloat);
315 }
316 if (sendtime) {
317 sprintf(pout, "time: %f ", coardstime); gaprnt(0, pout);
318 gawgdstime(outfile, &coardstime);
319 /* fwrite(&coardstime, sizeof(double), 1, outfile); */
320 }
321 gaprnt(0, "\n\t");
322
323 /* level independent data */
324 /* write data value and move ptr to next report simultaneously */
325 for (i = 0; i < numvars; i++) {
326 if (varlevels[i]) continue;
327 if (currpt[i] == NULL ||
328 currpt[i]->lat != ref->lat ||
329 currpt[i]->lon != ref->lon ||
330 currpt[i]->tim != ref->tim) {
331 gaprnt(0, "error: bad structure in result\n");
332 retval = 1;
333 goto cleanup;
334 }
335 if (sendind) {
336 sprintf(pout, "[%s: %f] ",
337 pcm->result[i].stn->pvar->abbrv, currpt[i]->val);
338 gaprnt(0, pout);
339 outFloat = currpt[i]->val;
340 gawgdsval(outfile, &outFloat);
341 }
342 currpt[i] = currpt[i]->rpt;
343 }
344
345 /* level dependent data */
346 /* write data for z-levels until we hit a new lat/lon/time */
347 if (levelstart < numvars) { /* if there are level-dep vars */
348 while (currpt[levelstart]&& /* and there are still more reports */
349 ref->lat == currpt[levelstart]->lat && /* and we haven't hit a */
350 ref->lon == currpt[levelstart]->lon && /* new lat/lon/time */
351 ref->tim == currpt[levelstart]->tim) {
352
353 levelref = currpt[levelstart];
354
355 gaprnt(0, "\n\tSOI ");
356 fwrite(startrec, sizeof(char), 4, outfile);
357
358 if (sendlev) {
359 sprintf(pout, "lev: %f ", levelref->lev); gaprnt(0, pout);
360 outFloat = levelref->lev;
361 gawgdsval(outfile, &outFloat);
362 }
363 /* write data value and move ptr to next report simultaneously */
364 for (i = levelstart; i < numvars; i++) {
365 if (!varlevels[i]) continue;
366 if (currpt[i] == NULL ||
367 currpt[i]->lat != ref->lat ||
368 currpt[i]->lon != ref->lon ||
369 currpt[i]->tim != ref->tim ||
370 currpt[i]->lev != levelref->lev) {
371 gaprnt(0, "error: bad structure in result\n");
372 retval = 1;
373 goto cleanup;
374 }
375 if (senddep) {
376 sprintf(pout, "[%s: %f] ", pcm->result[i].stn->pvar->abbrv,
377 currpt[i]->val); gaprnt(0, pout);
378 outFloat = currpt[i]->val;
379 gawgdsval(outfile, &outFloat);
380 }
381 currpt[i] = currpt[i]->rpt;
382 }
383 }
384 gaprnt(0, "\n\tEOS ");
385 fwrite(endrec, sizeof(char), 4, outfile);
386 }
387
388 gaprnt(0, "\n");
389 numreps++;
390 }
391
392 /* don't write the final EOS, so that time loops can be concatenated
393 * as a single sequence. */
394 /* gaprnt(0, "EOS\n"); */
395 /* fwrite(endrec, sizeof(char), 4, outfile); */
396
397 sprintf(pout, "WRITEGDS: %d reports x %d vars written as %d records\n",
398 pcm->result[0].stn->rnum, numvars, numreps); gaprnt(0, pout);
399
400 cleanup:
401 fclose(outfile);
402 free(varlevels);
403 free(currpt);
404 return;
405 }
406
407
gadprnt(struct gacmn * pcm)408 void gadprnt (struct gacmn *pcm) {
409 struct gastn *stn;
410 struct garpt *rpt;
411 struct gagrid *pgr;
412 float *gr;
413 int siz,i,j,k,lnum;
414
415 if (pcm->type[0] == 1) { /* Data type grid */
416 pgr = pcm->result[0].pgr;
417 siz = pgr->isiz*pgr->jsiz;
418 sprintf (pout,"Printing Grid -- %i Values -- Undef = %g\n",
419 siz, pgr->undef);
420 gaprnt(2,pout);
421 gr = pgr->grid;
422 lnum = 0;
423 for (i=0; i<siz; i++) {
424 if (pcm->prstr) {
425 if (*gr==pgr->undef && pcm->prudef) {
426 pout[0]='U'; pout[1]='n'; pout[2]='d';
427 pout[3]='e'; pout[4]='f'; pout[5]='\0';
428 } else {
429 sprintf (pout,pcm->prstr,*gr);
430 }
431 if (pcm->prbnum>0) {
432 j = 0;
433 while (pout[j]) j++;
434 for (k=0; k<pcm->prbnum; k++) {
435 pout[j] = ' ';
436 j++;
437 }
438 pout[j] = '\0';
439 }
440 gaprnt (2,pout);
441 } else {
442 sprintf (pout,"%g ",*gr);
443 gaprnt (2,pout);
444 }
445 lnum++;
446 if (lnum >= pcm->prlnum) {
447 gaprnt (2,"\n");
448 lnum = 0;
449 }
450 gr++;
451 }
452 if (lnum>0) gaprnt (2,"\n");
453 } else { /* Data type station */
454 stn = pcm->result[0].stn;
455 sprintf (pout,"Printing Stations -- %i Reports -- Undef = %g\n",
456 stn->rnum, stn->undef);
457 gaprnt(2,pout);
458 rpt = stn->rpt;
459 while (rpt) {
460 sprintf (pout,"%c%c%c%c%c%c%c%c %-8.6g %-8.6g %-8.6g \n",
461 rpt->stid[0], rpt->stid[1], rpt->stid[2], rpt->stid[3],
462 rpt->stid[4], rpt->stid[5], rpt->stid[6], rpt->stid[7],
463 rpt->lon,rpt->lat,rpt->lev);
464 gaprnt(2,pout);
465 if (pcm->prstr) {
466 sprintf(pout,pcm->prstr,rpt->val);
467 } else {
468 sprintf (pout,"%g ",rpt->val);
469 }
470 gaprnt(2,pout);
471 gaprnt(2,"\n");
472 rpt = rpt->rpt;
473 }
474 }
475 }
476
477 /* Write info and stats on data item to grads output stream */
478
gastts(struct gacmn * pcm)479 void gastts (struct gacmn *pcm) {
480 struct gastn *stn;
481 struct garpt *rpt;
482 struct gagrid *pgr;
483 struct dt dtim;
484 float (*conv) (float *, float);
485 float cint,cmin,cmax,rmin,rmax,*gr;
486 float sum,sumsqr,dum; /*mf mf*/
487 float gcntm1; /*mf mf*/
488 int i,j,ucnt,gcnt,gcnto,siz;
489 char lab[20];
490
491 if (pcm->type[0] == 1) { /* Data type grid */
492 pgr = pcm->result[0].pgr;
493 gaprnt(2,"Data Type = grid\n");
494 sprintf (pout,"Dimensions = %i %i\n",pgr->idim,pgr->jdim);
495 gaprnt(2,pout);
496 if (pgr->idim>-1) {
497 sprintf (pout,"I Dimension = %i to %i",pgr->dimmin[pgr->idim],
498 pgr->dimmax[pgr->idim]);
499 gaprnt(2,pout);
500
501 if (pgr->idim>-1 && pgr->ilinr==1) { /* Linear scaling info */
502 gaprnt(2," Linear");
503 if (pgr->idim==3) {
504 gr2t (pgr->ivals,pgr->dimmin[3],&dtim);
505 if (dtim.mn==0) gat2ch (&dtim,4,lab);
506 else gat2ch (&dtim,5,lab);
507 if (*(pgr->ivals+5)!=0) {
508 sprintf (pout," %s %gmo\n",lab,*(pgr->ivals+5));
509 } else {
510 sprintf (pout," %s %gmn\n",lab,*(pgr->ivals+6));
511 }
512 gaprnt (2,pout);
513 } else {
514 conv = pgr->igrab;
515 sprintf (pout," %g %g\n",conv(pgr->ivals,pgr->dimmin[pgr->idim]),
516 *(pgr->ivals));
517 gaprnt (2,pout);
518 }
519 }
520 if (pgr->idim>-1 && pgr->ilinr!=1) { /* Levels scaling info */
521 gaprnt(2," Levels");
522 conv = pgr->igrab;
523 for (i=pgr->dimmin[pgr->idim]; i<=pgr->dimmax[pgr->idim]; i++) {
524 sprintf (pout," %g",conv(pgr->ivals,i));
525 gaprnt (2,pout);
526 }
527 gaprnt (2,"\n");
528 }
529 } else {
530 gaprnt(2,"I Dimension = -999 to -999\n");
531 }
532 if (pgr->jdim>-1) {
533 sprintf (pout,"J Dimension = %i to %i",pgr->dimmin[pgr->jdim],
534 pgr->dimmax[pgr->jdim]);
535 gaprnt(2,pout);
536
537 if (pgr->jdim>-1 && pgr->jlinr==1) { /* Linear scaling info */
538 gaprnt(2," Linear");
539 if (pgr->jdim==3) {
540 gr2t (pgr->jvals,pgr->dimmin[3],&dtim);
541 if (dtim.mn==0) gat2ch (&dtim,4,lab);
542 else gat2ch (&dtim,5,lab);
543 if (*(pgr->jvals+5)!=0) {
544 sprintf (pout," %s %gmo\n",lab,*(pgr->jvals+5));
545 } else {
546 sprintf (pout," %s %gmn\n",lab,*(pgr->jvals+6));
547 }
548 gaprnt (2,pout);
549 } else {
550 conv = pgr->jgrab;
551 sprintf (pout," %g %g\n",conv(pgr->jvals,pgr->dimmin[pgr->jdim]),
552 *(pgr->jvals));
553 gaprnt (2,pout);
554 }
555 }
556 if (pgr->jdim>-1 && pgr->jlinr!=1) { /* Levels scaling info */
557 gaprnt(2," Levels");
558 conv = pgr->jgrab;
559 for (i=pgr->dimmin[pgr->jdim]; i<=pgr->dimmax[pgr->jdim]; i++) {
560 sprintf (pout," %g",conv(pgr->jvals,i));
561 gaprnt (2,pout);
562 }
563 gaprnt (2,"\n");
564 }
565 } else {
566 gaprnt(2,"J Dimension = -999 to -999\n");
567 }
568 siz = pgr->isiz*pgr->jsiz;
569 sprintf (pout,"Sizes = %i %i %i\n",pgr->isiz,pgr->jsiz,siz);
570 gaprnt(2,pout);
571 sprintf (pout,"Undef value = %g\n",pgr->undef);
572 gaprnt(2,pout);
573 ucnt = 0; gcnt = 0; sum=0; sumsqr=0; /*mf mf*/
574 gr = pgr->grid;
575 for (i=0; i<siz; i++) {
576 if (*(gr+i)==pgr->undef) ucnt++;
577 else{
578 dum = *(gr+i);
579 sum += dum;
580 sumsqr += dum*dum;
581 gcnt++;
582 }
583 }
584 sprintf (pout,"Undef count = %i Valid count = %i\n",ucnt,gcnt);
585 gaprnt(2,pout);
586 if (pgr->idim>-1) {
587 gamnmx (pgr);
588 sprintf (pout,"Min, Max = %g %g\n",pgr->rmin,pgr->rmax);
589 gaprnt(2,pout);
590 cint = 0.0;
591 gacsel (pgr->rmin,pgr->rmax,&cint,&cmin,&cmax);
592
593 if (pgr->jdim==-1) {
594 cmin = cmin - cint*2.0;
595 cmax = cmax + cint*2.0;
596 }
597 if (cint==0.0 || cint==pgr->undef) {
598 cmin = pgr->rmin-5.0;
599 cmax = pgr->rmax+5.0;
600 cint = 1.0;
601 }
602 sprintf (pout,"Cmin, cmax, cint = %g %g %g\n",cmin,cmax,cint);
603 gaprnt(2,pout);
604 /*mf------------mf*/
605 gcntm1=gcnt-1;
606 if(gcntm1<=0) gcntm1=1;
607 gcnto=gcnt;
608 if(gcnt<=0) gcnt=1;
609 sprintf (pout,"Stats[sum,sumsqr,root(sumsqr),n]: %g %g %g %d\n",sum,sumsqr,sqrt(sumsqr),gcnto);
610 gaprnt(2,pout);
611 sprintf (pout,"Stats[(sum,sumsqr,root(sumsqr))/n]: %g %g %g\n",sum/gcnt,sumsqr/gcnt,sqrt(sumsqr/gcnt));
612 gaprnt(2,pout);
613 sprintf (pout,"Stats[(sum,sumsqr,root(sumsqr))/(n-1)]: %g %g %g\n",sum/gcntm1,sumsqr/gcntm1,sqrt(sumsqr/gcntm1));
614 gaprnt(2,pout);
615 dum=(sumsqr/gcnt)-((sum/gcnt)*(sum/gcnt));
616 if(dum>0){
617 sprintf (pout,"Stats[(sigma,var)(n)]: %g %g\n",sqrt(dum),dum);
618 } else {
619 sprintf (pout,"Stats[(sigma,var)(n)]: %g %g\n",0.0,0.0);
620 }
621 gaprnt(2,pout);
622 dum=dum*(gcnt/gcntm1);
623 if(dum>0) {
624 sprintf (pout,"Stats[(sigma,var)(n-1)]: %g %g\n",sqrt(dum),dum);
625 } else {
626 sprintf (pout,"Stats[(sigma,var)(n-1)]: %g %g\n",0.0,0.0);
627 }
628 gaprnt(2,pout);
629 /*mf------------mf*/
630 } else {
631 sprintf (pout,"Min, Max = %g %g\n",pgr->rmin,pgr->rmin);
632 gaprnt(2,pout);
633 }
634
635 } else { /* Data type station */
636 gaprnt(2,"Data Type = station\n");
637 stn = pcm->result[0].stn;
638 sprintf (pout,"Dimensions = %i %i\n",stn->idim,stn->jdim);
639 gaprnt(2,pout);
640 if (stn->idim>-1) {
641 if (stn->idim!=3) {
642 sprintf (pout,"I Dimension = %g to %g\n",stn->dmin[stn->idim],stn->dmax[stn->idim]);
643 gaprnt(2,pout);
644 } else {
645 sprintf (pout,"I Dimension = %i to %i\n",stn->tmin, stn->tmax);
646 gaprnt(2,pout);
647 }
648 } else {
649 gaprnt(2,"I Dimension = -999 to -999\n");
650 }
651 if (stn->jdim>-1) {
652 if (stn->jdim!=3) {
653 sprintf (pout,"J Dimension = %g to %g\n",stn->dmin[stn->jdim],stn->dmax[stn->jdim]);
654 gaprnt(2,pout);
655 } else {
656 sprintf (pout,"J Dimension = %i to %i\n",stn->tmin, stn->tmax);
657 gaprnt(2,pout);
658 }
659 } else {
660 gaprnt(2,"J Dimension = -999 to -999\n");
661 }
662 sprintf (pout,"Stn count = %i\n",stn->rnum);
663 gaprnt(2,pout);
664 sprintf (pout,"Undef value = %g\n",stn->undef);
665 gaprnt(2,pout);
666 ucnt = 0; gcnt = 0; sum=0; sumsqr=0; /*mf mf*/
667 rmin = 9e33;
668 rmax = -9e33;
669 rpt = stn->rpt;
670
671 while (rpt) {
672 if (rpt->val==stn->undef) ucnt++;
673 else {
674 gcnt++;
675 dum = rpt->val;
676 sum += dum;
677 sumsqr += dum*dum;
678 if (rpt->val<rmin) rmin = rpt->val;
679 if (rpt->val>rmax) rmax = rpt->val;
680 }
681 rpt = rpt->rpt;
682 }
683 gcntm1=gcnt-1;
684 if(gcntm1 <=0) gcntm1=1;
685 if (gcnt==0) {
686 rmin = stn->undef;
687 rmax = stn->undef;
688 }
689
690 sprintf (pout,"Undef count = %i Valid count = %i \n",
691 ucnt,gcnt);
692 gaprnt(2,pout);
693 sprintf (pout,"Min, Max = %g %g\n",rmin,rmax);
694 gaprnt(2,pout);
695 cint = 0.0;
696
697 gacsel (rmin,rmax,&cint,&cmin,&cmax);
698 if (stn->jdim==-1) {
699 cmin = cmin - cint*2.0;
700 cmax = cmax + cint*2.0;
701 }
702 if (cint==0.0 || cint==stn->undef) {
703 cmin = rmin-5.0;
704 cmax = rmax+5.0;
705 cint = 1.0;
706 }
707 sprintf (pout,"Cmin, cmax, cint = %g %g %g\n",cmin,cmax,cint);
708 gaprnt(2,pout);
709
710 /*mf------------mf*/
711 gcntm1=gcnt-1;
712 if(gcntm1<=0) gcntm1=1;
713 gcnto=gcnt;
714 if(gcnt<=0) gcnt=1;
715 sprintf (pout,"Stats[sum,sumsqr,root(sumsqr),n]: %g %g %g %d\n",sum,sumsqr,sqrt(sumsqr),gcnto);
716 gaprnt(2,pout);
717 sprintf (pout,"Stats[(sum,sumsqr,root(sumsqr))/n)]: %g %g %g\n",sum/gcnt,sumsqr/gcnt,sqrt(sumsqr/gcnt));
718 gaprnt(2,pout);
719 sprintf (pout,"Stats[(sum,sumsqr,root(sumsqr))/(n-1))]: %g %g %g\n",sum/gcntm1,sumsqr/gcntm1,sqrt(sumsqr/gcntm1));
720 gaprnt(2,pout);
721 dum=(sumsqr/gcnt)-((sum/gcnt)*(sum/gcnt));
722 if(dum>0){
723 sprintf (pout,"Stats[(sigma,var)(n)]: %g %g\n",sqrt(dum),dum);
724 } else {
725 sprintf (pout,"Stats[(sigma,var)(n)]: %g %g\n",0.0,0.0);
726 }
727 gaprnt(2,pout);
728 dum=dum*(gcnt/gcntm1);
729 if(dum>0) {
730 sprintf (pout,"Stats[(sigma,var)(n-1)]: %g %g\n",sqrt(dum),dum);
731 } else {
732 sprintf (pout,"Stats[(sigma,var)(n-1)]: %g %g\n",0.0,0.0);
733 }
734 gaprnt(2,pout);
735
736 if(pcm->stnprintflg) {
737 sprintf (pout,"Printing station values: #obs = %d\n",gcnt);
738 gaprnt(2,pout);
739 gcnt=0;
740 ucnt=0;
741 rpt = stn->rpt;
742 sprintf (pout,"OB ID LON LAT LEV VAL\n");
743 gaprnt(2,pout);
744 while (rpt) {
745 if (rpt->val==stn->undef) ucnt++;
746 else {
747 gcnt++;
748 sprintf (pout,"%-5i %.8s %-8.6g %-8.6g %-8.6g %-8.6g\n",
749 gcnt,rpt->stid,rpt->lon,rpt->lat,rpt->lev,rpt->val);
750 gaprnt(2,pout);
751 }
752 rpt = rpt->rpt;
753 }
754 }
755
756 /*mf------------mf*/
757
758 }
759 }
760
761 /* Special routine -- find closest station to an X,Y position. */
762
gafstn(struct gacmn * pcm)763 void gafstn (struct gacmn *pcm) {
764 struct gastn *stn;
765 struct garpt *rpt, *srpt;
766 struct gagrid *pgr;
767 float x,y,xpos,ypos,r,d,rlon;
768
769 if (pcm->numgrd<3 || pcm->type[0]!=0 || pcm->type[1]!=1 ||
770 pcm->type[2]!=1) {
771 gaprnt (0,"Error: Invalid data types for findstn\n");
772 return;
773 }
774
775 gamscl (pcm); /* Do map level scaling */
776 stn = pcm->result[0].stn;
777 pgr = pcm->result[1].pgr;
778 xpos = pgr->rmin;
779 pgr = pcm->result[2].pgr;
780 ypos = pgr->rmin;
781 rpt = stn->rpt;
782 srpt = NULL;
783 r = 1.0e30;
784 while (rpt) {
785 rlon = rpt->lon;
786 if (rlon<pcm->dmin[0]) rlon+=360.0;
787 if (rlon>pcm->dmax[0]) rlon-=360.0;
788 gxconv (rlon,rpt->lat,&x,&y,2);
789 d = hypot(x-xpos,y-ypos);
790 if (d<r) {
791 r = d;
792 srpt = rpt;
793 }
794 rpt = rpt->rpt;
795 }
796 if (srpt) {
797 srpt->stid[7] = '\0';
798 sprintf (pout,"%s %g %g %g\n",srpt->stid,srpt->lon,srpt->lat,r);
799 gaprnt(2,pout);
800 } else gaprnt (2,"No stations found\n");
801 gagsav (21,pcm,NULL);
802 }
803
804 /* Plot weather symbols at station locations. */
805
gawsym(struct gacmn * pcm)806 void gawsym (struct gacmn *pcm) {
807 struct gastn *stn;
808 struct garpt *rpt;
809 float rlon,x,y,x1,x2,y1,y2,scl;
810 int i;
811
812 gamscl (pcm); /* Do map level scaling */
813 gawmap (pcm, 1); /* Draw map */
814 pcm->xdim = 0; pcm->ydim = 1;
815 gafram (pcm);
816 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
817
818 stn = pcm->result[0].stn;
819 rpt = stn->rpt;
820 gxwide (pcm->cthick);
821 if (pcm->ccolor<0) gxcolr(1);
822 else gxcolr (pcm->ccolor);
823 while (rpt!=NULL) {
824 if (rpt->val != stn->undef) {
825 rlon = rpt->lon;
826 if (rlon<pcm->dmin[0]) rlon+=360.0;
827 if (rlon>pcm->dmax[0]) rlon-=360.0;
828 if (rlon>pcm->dmin[0]&&rlon<pcm->dmax[0]&&
829 rpt->lat>pcm->dmin[1]&&rpt->lat<pcm->dmax[1]) {
830 i = (int)(rpt->val+0.1);
831 if (i>0&&i<42) {
832 gxconv (rlon,rpt->lat,&x,&y,2);
833 scl = pcm->digsiz*1.5;
834 /*
835 gxcolr(0);
836 x1 = x - (sxwid[i-1]*0.7*scl);
837 x2 = x + (sxwid[i-1]*0.7*scl);
838 y1 = y + (symin[i-1]*scl*1.3);
839 y2 = y + (symax[i-1]*scl*1.2);
840 gxrecf (x1,x2,y1,y2);
841 */
842 gxwide (pcm->cthick);
843 if (pcm->wxopt==1) {
844 wxsym (i, x, y, scl, pcm->ccolor, pcm->wxcols);
845 } else {
846 gxmark (i,x,y,scl);
847 }
848 }
849 }
850 }
851 rpt = rpt->rpt;
852 }
853 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
854 gxcolr(pcm->anncol);
855 gxwide(pcm->annthk-3);
856 if (pcm->pass==0 && pcm->grdsflg)
857 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
858 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
859 gxstyl(1);
860 gaaxpl(pcm,0,1);
861 gagsav(12,pcm,NULL);
862 }
863
864 /* Plot colorized markers at station locations */
865
gasmrk(struct gacmn * pcm)866 void gasmrk (struct gacmn *pcm) {
867 struct gastn *stn;
868 struct garpt *rpt;
869 float rlon,x,y;
870 int i,icnst;
871 int cnt; /* qqq */
872 char lab[20];
873 float cwid,sizstid;
874 int len;
875
876 gamscl (pcm); /* Do map level scaling */
877 gawmap (pcm, 1); /* Draw map */
878 pcm->xdim = 0; pcm->ydim = 1;
879 gafram (pcm);
880 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
881
882 stn = pcm->result[0].stn;
883 gasmnmx (stn);
884
885 /* qqq
886 printf ("qqq smin, smax = %g %g \n",stn->smin,stn->smax);
887 printf("qqq cmark %i %i %d\n",pcm->cmark,pcm->pass,pcm->pfid->type);
888 printf("qqq smin smax %f %f %i\n",stn->smin,stn->smax,pcm->ccolor);
889 qqq */
890
891
892 if (stn->smin == stn->undef || stn->smax == stn->undef) return;
893 gaselc (pcm,stn->smin,stn->smax);
894 rpt = stn->rpt;
895 gxwide (pcm->cthick);
896
897 icnst=0;
898 if(stn->smin == stn->smax) icnst=1;
899
900 /* qqq
901 printf("qqq icnst %i %i\n",icnst,pcm->ccolor);
902 qqq */
903
904 if (pcm->ccolor<0 && icnst) {
905 pcm->ccolor=1;
906 gxcolr(pcm->ccolor);
907 } else if (pcm->ccolor<0) {
908 gxcolr(1);
909 } else {
910 gxcolr (pcm->ccolor);
911 }
912
913 cnt = 0; /* qqq */
914 sizstid=pcm->digsiz*0.65;
915
916 /*mf 981214 - default cmark is set to closed circle in gauser.c */
917
918 while (rpt!=NULL) {
919 cnt++;
920 if (rpt->val != stn->undef) {
921
922 /* qqq
923 if (rpt->val < -300.0) {
924 printf ("qqq %g %i %c%c%c%c%c\n",rpt->val,cnt,
925 rpt->stid[0],rpt->stid[1],rpt->stid[2],rpt->stid[3],rpt->stid[4]);
926 }
927 qqq */
928
929 rlon = rpt->lon;
930 if (rlon<pcm->dmin[0]) rlon+=360.0;
931 if (rlon>pcm->dmax[0]) rlon-=360.0;
932 if (rlon>pcm->dmin[0]&&rlon<pcm->dmax[0]&&
933 rpt->lat>pcm->dmin[1]&&rpt->lat<pcm->dmax[1]) {
934 gxconv (rlon,rpt->lat,&x,&y,2);
935 i = gashdc (pcm,rpt->val);
936
937 /*mf 981214 test if constant grid, if so use user ccolor */
938
939 if (icnst) {
940 gxcolr(pcm->ccolor);
941 } else {
942 gxcolr (i);
943 }
944
945 gxmark (pcm->cmark,x,y,pcm->digsiz*0.5);
946
947 /* 20000723 - stn id plot */
948
949 if (pcm->stidflg) {
950 gxmark (1,x,y,pcm->digsiz*0.5);
951 getwrd (lab,rpt->stid,8);
952 len = strlen(lab);
953 cwid = 0.1;
954 gxchln (lab,len,sizstid,&cwid);
955 x = x-cwid*0.5;
956 y = y-(sizstid*1.7);
957 if (pcm->ccolor!=0) {
958 gxcolr (gxqbck());
959 gxrecf (x-0.01,x+cwid+0.01,y-0.01,y+sizstid+0.01);
960 }
961 if (pcm->ccolor<0) gxcolr(1);
962 else gxcolr (pcm->ccolor);
963 gxchpl (lab,len,x,y,sizstid,sizstid,0.0);
964 }
965
966 }
967 }
968 rpt = rpt->rpt;
969 }
970 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
971 gxcolr(pcm->anncol);
972 gxwide(pcm->annthk-3);
973 if (pcm->pass==0 && pcm->grdsflg)
974 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
975 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
976 gxstyl(1);
977 gaaxpl(pcm,0,1);
978 gagsav(12,pcm,NULL);
979 }
980
981 /* Plot station values as either one number centered on the station
982 location or as two numbers above and below a crosshair, or as
983 a wind barb if specified by user */
984
gapstn(struct gacmn * pcm)985 void gapstn (struct gacmn *pcm) {
986 struct gastn *stn, *stn2;
987 int i,len,flag,hemflg;
988 struct garpt *rpt, *rpt2;
989 float x,y,rlon,dir,spd,umax,vmax,vscal,cwid;
990 char lab[20];
991
992 gamscl (pcm); /* Do map level scaling */
993 gawmap (pcm, 1); /* Draw map */
994 pcm->xdim = 0;
995 pcm->ydim = 1;
996 gafram(pcm);
997 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
998
999 /* Plot barb or vector at station location */
1000
1001 if (pcm->numgrd>1 && (pcm->goutstn==2 || pcm->goutstn==6)) {
1002 stn = pcm->result[0].stn;
1003 stn2 = pcm->result[1].stn;
1004 if (pcm->goutstn==6) { /* Get vector scaling */
1005 if (!pcm->arrflg) {
1006 rpt = stn->rpt;
1007 umax = -9.99e33;
1008 while (rpt) {
1009 if (umax<fabs(rpt->val)) umax = fabs(rpt->val);
1010 rpt = rpt->rpt;
1011 }
1012 rpt = stn2->rpt;
1013 vmax = -9.99e33;
1014 while (rpt) {
1015 if (vmax<fabs(rpt->val)) vmax = fabs(rpt->val);
1016 rpt = rpt->rpt;
1017 }
1018 vscal = hypot(umax,vmax);
1019 x = floor(log10(vscal));
1020 y = floor(vscal/pow(10.0,x));
1021 vscal = y * pow(10.0,x);
1022 pcm->arrsiz = 0.5;
1023 pcm->arrmag = vscal;
1024 } else {
1025 vscal = pcm->arrmag;
1026 }
1027 pcm->arrflg = 1;
1028 }
1029 rpt = stn->rpt;
1030 if (pcm->ccolor<0) gxcolr(1);
1031 else gxcolr (pcm->ccolor);
1032 gxwide (pcm->cthick);
1033 while (rpt!=NULL) {
1034 if (rpt->val != stn->undef) {
1035 rpt2 = stn2->rpt;
1036 while (rpt2!=NULL) {
1037 if (rpt2->val != stn->undef && rpt->lat==rpt2->lat &&
1038 rpt->lon==rpt2->lon) {
1039 rlon = rpt->lon;
1040 if (rlon<pcm->dmin[0]) rlon+=360.0;
1041 if (rlon>pcm->dmax[0]) rlon-=360.0;
1042 if (rlon>pcm->dmin[0]&&rlon<pcm->dmax[0]&&
1043 rpt->lat>pcm->dmin[1]&&rpt->lat<pcm->dmax[1]) {
1044 gxconv (rlon,rpt->lat,&x,&y,2);
1045 if (rpt2->val==0.0 && rpt->val==0.0) dir = 0.0;
1046 else dir = atan2(rpt2->val,rpt->val) +
1047 gxaarw(rpt->lon,rpt->lat);
1048 spd = hypot(rpt->val,rpt2->val);
1049 if (pcm->goutstn==2) {
1050 hemflg = 0;
1051 if (pcm->hemflg == 1) hemflg = 1;
1052 else if (pcm->hemflg == 0) hemflg = 0;
1053 else if (rpt->lat<0.0) hemflg = 1;
1054 gabarb (x, y, pcm->digsiz*3.5, pcm->digsiz*2.0,
1055 pcm->digsiz*0.25, dir, spd, hemflg);
1056 gxmark (2,x,y,pcm->digsiz*0.5);
1057 } else {
1058 if (vscal>0.0) {
1059 gaarrw (x, y, dir, pcm->arrsiz*spd/vscal, pcm->ahdsiz);
1060 } else {
1061 gaarrw (x, y, dir, pcm->arrsiz, pcm->ahdsiz);
1062 }
1063 }
1064 }
1065 break;
1066 }
1067 rpt2 = rpt2->rpt;
1068 }
1069 }
1070 rpt = rpt->rpt;
1071 }
1072
1073 /* Plot number at station location */
1074
1075 } else {
1076 gxwide (pcm->cthick);
1077 stn = pcm->result[0].stn;
1078 rpt = stn->rpt;
1079 flag=0;
1080 if (pcm->numgrd>1 || pcm->stidflg) flag = 1;
1081 while (rpt!=NULL) {
1082 if (rpt->val != stn->undef) {
1083 rlon = rpt->lon;
1084 if (rlon<pcm->dmin[0]) rlon+=360.0;
1085 if (rlon>pcm->dmax[0]) rlon-=360.0;
1086 if (rlon>pcm->dmin[0]&&rlon<pcm->dmax[0]&&
1087 rpt->lat>pcm->dmin[1]&&rpt->lat<pcm->dmax[1]) {
1088 gxconv (rlon,rpt->lat,&x,&y,2);
1089 if (flag) gxmark (1,x,y,pcm->digsiz*0.5);
1090 sprintf (lab,"%.*f",pcm->dignum,rpt->val);
1091 len = strlen(lab);
1092 cwid = 0.1;
1093 gxchln (lab,len,pcm->digsiz,&cwid);
1094 x = x-cwid*0.5;
1095 if (flag) {
1096 y = y+(pcm->digsiz*0.7);
1097 } else {
1098 y = y-(pcm->digsiz*0.5);
1099 }
1100 if (pcm->ccolor!=0) {
1101 gxcolr (gxqbck());
1102 gxrecf (x-0.01,x+cwid+0.01,y-0.01,y+pcm->digsiz+0.01);
1103 }
1104 if (pcm->ccolor<0) gxcolr(1);
1105 else gxcolr (pcm->ccolor);
1106 gxchpl (lab,len,x,y,pcm->digsiz,pcm->digsiz,0.0);
1107 }
1108 }
1109 rpt=rpt->rpt;
1110 }
1111 if ( (flag && pcm->type[1]==0) || pcm->stidflg) {
1112 if (!pcm->stidflg) stn = pcm->result[1].stn;
1113 rpt = stn->rpt;
1114 while (rpt!=NULL) {
1115 rlon = rpt->lon;
1116 if (rlon<pcm->dmin[0]) rlon+=360.0;
1117 if (rlon>pcm->dmax[0]) rlon-=360.0;
1118 if (rlon>pcm->dmin[0]&&rlon<pcm->dmax[0]&&
1119 rpt->lat>pcm->dmin[1]&&rpt->lat<pcm->dmax[1]) {
1120 gxconv (rlon,rpt->lat,&x,&y,2);
1121 gxmark (1,x,y,pcm->digsiz*0.5);
1122 if (pcm->stidflg) {
1123 getwrd (lab,rpt->stid,8);
1124 } else {
1125 sprintf (lab,"%.*f",pcm->dignum,rpt->val);
1126 }
1127 len = strlen(lab);
1128 cwid = 0.1;
1129 gxchln (lab,len,pcm->digsiz,&cwid);
1130 x = x-cwid*0.5;
1131 y = y-(pcm->digsiz*1.7);
1132 if (pcm->ccolor!=0) {
1133 gxcolr (gxqbck());
1134 gxrecf (x-0.01,x+cwid+0.01,y-0.01,y+pcm->digsiz+0.01);
1135 }
1136 if (pcm->ccolor<0) gxcolr(1);
1137 else gxcolr (pcm->ccolor);
1138 gxchpl (lab,len,x,y,pcm->digsiz,pcm->digsiz,0.0);
1139 }
1140 rpt=rpt->rpt;
1141 }
1142 }
1143 }
1144 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
1145 gxcolr(pcm->anncol);
1146 gxwide(pcm->annthk-3);
1147 if (pcm->pass==0 && pcm->grdsflg)
1148 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
1149 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
1150 gxstyl(1);
1151 gaaxpl(pcm,0,1);
1152 gagsav (10,pcm,NULL);
1153 }
1154
1155 /* Draw wind barb at location x, y, direction dir, speed
1156 spd, with barb lengths of blen, pointer length of plen
1157 (before barbs added), starting at rad distance from the
1158 center point x,y. If calm, a circle is drawn at rad*1.2
1159 radius. Direction is direction wind blowing towards. */
1160
gabarb(float x,float y,float plen,float blen,float rad,float dir,float spd,int hemflg)1161 void gabarb (float x, float y, float plen, float blen, float rad,
1162 float dir, float spd, int hemflg) {
1163 float bgap,padd,var,a70,cosd70,sind70;
1164 float xp1,yp1,xp2,yp2,xp3,yp3;
1165 int flag,flg2;
1166
1167 dir = dir + 3.14159; /* Want direction wind blowing from */
1168
1169 /* Calculate added length due to barbs */
1170
1171 bgap = blen*0.3;
1172 padd = 0.0;
1173 var = spd;
1174 if (var<0.01) {
1175 rad*=2.0;
1176 if (rad+blen*0.3>rad*1.4) gxmark(2,x,y,rad+blen*0.3);
1177 else gxmark (2,x,y,rad*1.4);
1178 return;
1179 } else {
1180 if (var<2.5) padd = bgap;
1181 while (var>=47.5) { var-=50.0; padd+=bgap*1.6;}
1182 while (var>=7.5) { var-=10.0; padd+=bgap;}
1183 if (var>=2.5) padd+=bgap;
1184 }
1185 plen+=padd;
1186
1187 /* Draw pointer */
1188
1189 xp1 = x + plen*cos(dir);
1190 yp1 = y + plen*sin(dir);
1191 xp2 = x + rad*cos(dir);
1192 yp2 = y + rad*sin(dir);
1193 gxplot (xp2,yp2,3);
1194 gxplot (xp1,yp1,2);
1195
1196 /* Start out at the end of the pointer and add barbs
1197 til we run out of barbs to add. */
1198
1199 var = spd;
1200 a70 = 70.0*3.14159/180.0;
1201 if (hemflg) {
1202 cosd70 = cos(dir+a70);
1203 sind70 = sin(dir+a70);
1204 } else {
1205 cosd70 = cos(dir-a70);
1206 sind70 = sin(dir-a70);
1207 }
1208 flag = 1;
1209 flg2 = 0;
1210 while (var>=47.5) {
1211 xp1 = x + plen*cos(dir);
1212 yp1 = y + plen*sin(dir);
1213 xp2 = xp1 + blen*cosd70;
1214 yp2 = yp1 + blen*sind70;
1215 xp3 = x + (plen-bgap*1.45)*cos(dir);
1216 yp3 = y + (plen-bgap*1.45)*sin(dir);
1217 gxplot (xp1,yp1,3);
1218 gxplot (xp2,yp2,2);
1219 gxplot (xp3,yp3,2);
1220 plen = plen - bgap*1.6;
1221 var-=50.0;
1222 flg2 = 1;
1223 }
1224 while (var>=7.5) {
1225 if (flg2) {plen-=bgap*0.7; flg2 = 0;}
1226 xp1 = x + plen*cos(dir);
1227 yp1 = y + plen*sin(dir);
1228 xp2 = xp1 + blen*cosd70;
1229 yp2 = yp1 + blen*sind70;
1230 gxplot (xp1,yp1,3);
1231 gxplot (xp2,yp2,2);
1232 plen-=bgap;
1233 var-=10.0;
1234 flag = 0;
1235 }
1236 if (var>=2.5) {
1237 if (flag) plen-=bgap;
1238 xp1 = x + plen*cos(dir);
1239 yp1 = y + plen*sin(dir);
1240 xp2 = xp1 + 0.5*blen*cosd70;
1241 yp2 = yp1 + 0.5*blen*sind70;
1242 gxplot (xp1,yp1,3);
1243 gxplot (xp2,yp2,2);
1244 }
1245 return;
1246 }
1247
1248 /* Plot station model using station data. Complexity of model
1249 depends on the number of result objects available.
1250 First two are assumed to be u and v winds, which are required
1251 for the station model to be plotted. */
1252
gapmdl(struct gacmn * pcm)1253 void gapmdl (struct gacmn *pcm) {
1254 struct gastn *stn, *stn1;
1255 struct garpt *rpt, *rpt1;
1256 struct gagrid *pgr;
1257 float vals[10],udefs[10];
1258 int i,num;
1259
1260 gamscl (pcm); /* Do map level scaling */
1261 gawmap (pcm, 1); /* Draw map */
1262 pcm->xdim = 0;
1263 pcm->ydim = 1;
1264 gafram (pcm);
1265 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
1266
1267 for (i=0; i<10; i++) {
1268 vals[i] = -999.9;
1269 udefs[i] = -999.9;
1270 }
1271
1272 num = pcm->numgrd;
1273 if (num>10) num=10;
1274
1275 /* Check validity of data */
1276
1277 for (i=0; i<pcm->numgrd; i++) {
1278 if (pcm->type[i]==1) {
1279 pgr = pcm->result[i].pgr;
1280 if (pgr->idim!=-1) {
1281 gaprnt (0,"Invalid data: Station data required\n");
1282 return;
1283 }
1284 udefs[i] = *(pgr->grid);
1285 } else {
1286 stn = pcm->result[i].stn;
1287 udefs[i] = stn->undef;
1288 }
1289 vals[i] = udefs[i];
1290 }
1291
1292 if (pcm->ccolor<0) gxcolr(1);
1293 else gxcolr (pcm->ccolor);
1294 gxwide (pcm->cthick);
1295
1296 /* Get info for each station */
1297
1298 stn1 = pcm->result[0].stn;
1299 rpt1 = stn1->rpt;
1300 while (rpt1) {
1301 vals[0] = rpt1->val;
1302 for (i=1; i<pcm->numgrd; i++) {
1303 if (pcm->type[i]==1) {
1304 pgr = pcm->result[i].pgr;
1305 vals[i] = *(pgr->grid);
1306 } else {
1307 stn = pcm->result[i].stn;
1308 rpt = stn->rpt;
1309 while (rpt) {
1310 if (rpt->lon==rpt1->lon && rpt->lat==rpt1->lat &&
1311 cmpch(rpt->stid,rpt1->stid,8)==0) break;
1312 rpt = rpt->rpt;
1313 }
1314 if (rpt) vals[i] = rpt->val;
1315 else vals[i] = udefs[i];
1316 }
1317 }
1318 gasmdl (pcm,rpt1,vals,udefs);
1319 rpt1 = rpt1->rpt;
1320 }
1321
1322 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
1323 gxcolr(pcm->anncol);
1324 gxwide(pcm->annthk-3);
1325 if (pcm->pass==0 && pcm->grdsflg)
1326 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
1327 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
1328 gxstyl(1);
1329 gaaxpl(pcm,0,1);
1330 gagsav (11,pcm,NULL);
1331 }
1332
1333 /* Plot an individual station model */
1334
gasmdl(struct gacmn * pcm,struct garpt * rpt,float * vals,float * udefs)1335 void gasmdl (struct gacmn *pcm, struct garpt *rpt,
1336 float *vals, float *udefs) {
1337 int num,icld,i,col,ivis,itop,ibot,ii,hemflg;
1338 float x,y,spd,dir,roff,scl,rad,t,digsiz,msize,plen,wrad;
1339 float xlo[7],xhi[7],ylo[7],yhi[7],orad,xv,yv;
1340 char ch[20],len;
1341
1342 col = pcm->ccolor;
1343 if (pcm->ccolor<0) col = 1;
1344 gxcolr(col);
1345
1346 num = pcm->numgrd;
1347
1348 /* If no winds, just return */
1349
1350 if (*vals==*udefs || *(vals+1)==*(udefs+1)) return;
1351
1352 /* Plot */
1353
1354 gxconv (rpt->lon,rpt->lat,&x,&y,2);
1355
1356 hemflg = 0;
1357 if (pcm->hemflg == 1) hemflg = 1;
1358 else if (pcm->hemflg == 0) hemflg = 0;
1359 else if (rpt->lat<0.0) hemflg = 1;
1360
1361 icld = (int)(floor(*(vals+6)+0.1));
1362 digsiz = pcm->digsiz;
1363 if (icld>19 && icld<26) msize = digsiz*1.4;
1364 else msize = digsiz;
1365 wrad = msize*0.5;
1366 if (*vals==0.0 && *(vals+1)==0.0) rad = msize*0.8;
1367 else rad = msize*0.65;
1368
1369
1370 /* Plot clouds: 20 = clr, 21 = sct, 22 = bkn, 23 = ovc, 24 = obsc,
1371 25 = missng */
1372
1373 icld = (int)(floor(*(vals+6)+0.1));
1374 if (col==0 && icld>19) {
1375 if (icld!=23) gxmark(3,x,y,msize);
1376 } else {
1377 if (icld==20) gxmark(2,x,y,msize);
1378 else if (icld==21) gxmark(10,x,y,msize);
1379 else if (icld==22) gxmark(11,x,y,msize);
1380 else if (icld==23) gxmark(3,x,y,msize);
1381 else if (icld==24) {
1382 gxmark(2,x,y,msize);
1383 roff = msize/2.8284;
1384 gxplot (x-roff,y+roff,3);
1385 gxplot (x+roff,y-roff,2);
1386 gxplot (x+roff,y+roff,3);
1387 gxplot (x-roff,y-roff,2);
1388 }
1389 else if (icld==25) {
1390 gxmark(2,x,y,msize);
1391 roff = msize*0.3;
1392 gxchpl ("`0M",3,x-roff*1.2,y-roff*0.9,roff*2.0,roff*2.0,0.0);
1393 }
1394 else {
1395 if (icld>0 && icld<9) gxmark(icld,x,y,msize);
1396 else gxmark(2,x,y,msize);
1397 }
1398 }
1399
1400 /* Plot wxsym */
1401
1402 gxwide (pcm->cthick+1);
1403 xlo[0] = -999.9;
1404 if (*(vals+7)!=*(udefs+7)) {
1405 i = (int)(*(vals+7)+0.1);
1406 if (i>0&&i<40) {
1407 scl = digsiz*3.0;
1408 xhi[0] = x - rad*1.1;
1409 xlo[0] = xhi[0] - (sxwid[i-1]*scl);
1410 ylo[0] = y + symin[i-1]*scl;
1411 yhi[0] = y + symax[i-1]*scl;
1412 wxsym (i,xlo[0]+(sxwid[i-1]*scl*0.5),y,scl,pcm->ccolor,pcm->wxcols);
1413 gxcolr(col);
1414 }
1415 }
1416
1417 /* Plot visibility */
1418
1419 gxwide(pcm->cthick);
1420 xlo[5] = -999.9;
1421 if (*(vals+8) != *(udefs+8)) {
1422 i = (int)((*(vals+8)*32.0)+0.5);
1423 ivis = i/32;
1424 itop = i - ivis*32;
1425 ibot = 32;
1426 while (itop!=0) {
1427 i = itop/2;
1428 if (i*2!=itop) break;
1429 itop = i;
1430 ibot = ibot/2;
1431 }
1432 yv = y-digsiz*0.5;
1433 if (xlo[0]<-990.0) xv = x - rad*1.6;
1434 else xv = xlo[0] - digsiz*0.2;
1435 xhi[5] = xv;
1436 ylo[5] = yv;
1437 yhi[5] = yv+digsiz;
1438 if (itop>0) {
1439 sprintf (ch,"%i",ibot);
1440 len = 0;
1441 while (*(ch+len)) len++;
1442 xv = xv - ((float)len)*digsiz*0.65;
1443 gxchpl (ch,len,xv,yv,digsiz*0.65,digsiz*0.65,0.0);
1444 sprintf (ch,"%i",itop);
1445 len = 0;
1446 while (*(ch+len)) len++;
1447 gxplot (xv-digsiz*0.4,yv,3);
1448 gxplot (xv+digsiz*0.1,yv+digsiz,2);
1449 xv = xv - ((float)len+0.4)*digsiz*0.65;
1450 gxchpl (ch,len,xv,yv+digsiz*0.35,digsiz*0.65,digsiz*0.65,0.0);
1451 }
1452 if (ivis>0 || (ivis==0 && itop==0)) {
1453 sprintf (ch,"%i",ivis);
1454 len = 0;
1455 while (*(ch+len)) len++;
1456 xv = xv - ((float)len)*digsiz;
1457 gxchpl (ch,len,xv,yv,digsiz,digsiz,0.0);
1458 }
1459 xlo[5] = xv;
1460 }
1461
1462 /* Plot temperature */
1463
1464 xlo[1] = -999.9;
1465 if (*(vals+2) != *(udefs+2)) {
1466 sprintf (ch,"%.*f",pcm->dignum,*(vals+2));
1467 len = 0;
1468 while (*(ch+len)) len++;
1469 if (xlo[0]>-999.0) ylo[1] = yhi[0]+digsiz*0.4;
1470 else if (xlo[5]>-999.0) ylo[1] = yhi[5]+digsiz*0.4;
1471 else ylo[1] = y + digsiz*0.3;
1472 t = ylo[1]-y;
1473 t = rad*rad - t*t;
1474 if (t<=0.0) xhi[1] = x-digsiz*0.3;
1475 else {
1476 t = sqrt(t);
1477 if (t<digsiz*0.3) t = digsiz*0.3;
1478 xhi[1] = x-t;
1479 }
1480 xlo[1] = xhi[1] - (float)len * digsiz;
1481 yhi[1] = ylo[1] + digsiz;
1482 gxchpl (ch,len,xlo[1],ylo[1],digsiz,digsiz,0.0);
1483 }
1484
1485 /* Plot dewpoint */
1486
1487 xlo[2] = -999.9;
1488 if (*(vals+3) != *(udefs+3)) {
1489 sprintf (ch,"%.*f",pcm->dignum,*(vals+3));
1490 len = 0;
1491 while (*(ch+len)) len++;
1492 if (xlo[0]>-999.0) yhi[2] = ylo[0]-digsiz*0.4;
1493 else if (xlo[5]>-999.0) yhi[2] = ylo[5]-digsiz*0.4;
1494 else yhi[2] = y - digsiz*0.3;
1495 t = y - yhi[2];
1496 t = rad*rad - t*t;
1497 if (t<=0.0) xhi[2] = x-digsiz*0.3;
1498 else {
1499 t = sqrt(t);
1500 if (t<digsiz*0.3) t = digsiz*0.3;
1501 xhi[2] = x-t;
1502 }
1503 xlo[2] = xhi[2] - (float)len * digsiz;
1504 ylo[2] = yhi[2] - digsiz;
1505 gxchpl (ch,len,xlo[2],ylo[2],digsiz,digsiz,0.0);
1506 }
1507
1508 /* Plot Pressure */
1509
1510 xlo[3] = -999.9;
1511 if (*(vals+4) != *(udefs+4)) {
1512 i = (int)(*(vals+4)+0.5);
1513 ii = 0;
1514 if (pcm->mdldig3) {
1515 sprintf (ch,"%i",i+10000);
1516 len = 0;
1517 while (*(ch+len)) len++;
1518 ii = len-3;
1519 len = 3;
1520 } else {
1521 sprintf (ch,"%.*f",pcm->dignum,*(vals+4));
1522 len = 0;
1523 while (*(ch+len)) len++;
1524 }
1525 xlo[3] = x + rad*0.87;
1526 ylo[3] = y + rad*0.5;
1527 xhi[3] = xlo[3] + digsiz*(float)len;
1528 yhi[3] = ylo[3] + digsiz;
1529 gxchpl (ch+ii,len,xlo[3],ylo[3],digsiz,digsiz,0.0);
1530 }
1531
1532 /* Plot pressure tendency */
1533
1534 xlo[4] = -999.9;
1535 if (*(vals+5) != *(udefs+5)) {
1536 if (*(vals+5)>0.0) sprintf (ch,"+%.*f",pcm->dignum,*(vals+5));
1537 else sprintf (ch,"%.*f",pcm->dignum,*(vals+5));
1538 len = 0;
1539 while (*(ch+len)) len++;
1540 xlo[4] = x + rad;
1541 if (xlo[3]>-990.0) yhi[4] = ylo[3]-digsiz*0.4;
1542 else yhi[4] = y + digsiz*0.5;
1543 xhi[4] = xlo[4] + digsiz*(float)len;
1544 ylo[4] = yhi[4] - digsiz;
1545 gxchpl (ch,len,xlo[4],ylo[4],digsiz,digsiz,0.0);
1546 }
1547
1548 /* xxxx plot stid lower right */
1549
1550 xlo[6] = -999.9;
1551 if (pcm->stidflg) {
1552 len = 4;
1553 if (xlo[4]>-990.0) yhi[6] = ylo[4]-digsiz*0.4;
1554 else yhi[6] = y - rad;
1555 if (xlo[2]>-990.0) xlo[6] = xhi[2]+0.6*digsiz;
1556 else xlo[6] = x - 1.2*digsiz;
1557 xhi[6] = xlo[6] + 0.6*digsiz*(float)len;
1558 ylo[6] = yhi[6] - 0.6*digsiz;
1559 gxchpl (rpt->stid,len,xlo[6],ylo[6],0.6*digsiz,0.6*digsiz,0.0);
1560 }
1561
1562 /* Plot wind barb */
1563
1564 if (*vals==0.0 && *(vals+1)==0.0) dir = 0.0;
1565 else dir = atan2(*(vals+1),*vals) + gxaarw(rpt->lon,rpt->lat);
1566 spd = hypot(*vals,*(vals+1));
1567
1568 orad = wndexit (spd*cos(dir),spd*sin(dir), x, y, rad, xlo, xhi, ylo, yhi);
1569 if (orad<-990.0) {
1570 plen = pcm->digsiz*3.5;
1571 } else {
1572 plen = pcm->digsiz*0.5+orad;
1573 if (plen<pcm->digsiz*3.5) plen = pcm->digsiz*3.5;
1574 if (plen>pcm->digsiz*6.0) plen = orad;
1575 wrad = orad;
1576 }
1577 gabarb (x, y, plen, pcm->digsiz*2.5, wrad, dir, spd, hemflg);
1578 }
1579
1580 /* Find exit radius for the wind barb */
1581
wndexit(float uu,float vv,float x,float y,float rad,float * xlo,float * xhi,float * ylo,float * yhi)1582 float wndexit (float uu, float vv, float x, float y, float rad,
1583 float *xlo, float *xhi, float *ylo, float *yhi) {
1584 float u,v,xn,yn,orad,fuzz,fuzz2;
1585
1586 u = -1.0*uu;
1587 v = -1.0*vv;
1588 orad = -999.9;
1589 fuzz = rad*0.25;
1590 fuzz2 = fuzz*0.5;
1591
1592 if (u==0.0 && v==0.0) return(orad);
1593 if (u<0.0) {
1594 if (v>0.0 && xlo[1]>-990.0) {
1595 yhi[1] = yhi[1]+fuzz;
1596 xn = x + (yhi[1]-y)*u/v;
1597 if (xn>xlo[1]-fuzz && xn<xhi[1]+fuzz2) {
1598 orad = hypot(xn-x, yhi[1]-y);
1599 return (orad);
1600 }
1601 xlo[1] = xlo[1] - fuzz2;
1602 yn = y + (xlo[1]-x)*v/u;
1603 if (yn>ylo[1]-fuzz && yn<yhi[1]+fuzz) {
1604 orad = hypot(xlo[1]-x,yn-y);
1605 return (orad);
1606 }
1607 }
1608 if (v<0.0 && xlo[2]>-990.0) {
1609 ylo[2] = ylo[2]-fuzz;
1610 xn = x + (ylo[2]-y)*u/v;
1611 if (xn>xlo[2]-fuzz && xn<xhi[2]+fuzz2) {
1612 orad = hypot(xn-x, ylo[2]-y);
1613 return (orad);
1614 }
1615 xlo[2] = xlo[2] - fuzz2;
1616 yn = y + (xlo[2]-x)*v/u;
1617 if (yn>ylo[2]-fuzz && yn<yhi[2]+fuzz) {
1618 orad = hypot(xlo[2]-x,yn-y);
1619 return (orad);
1620 }
1621 }
1622 if (xlo[5]>-990.0) {
1623 xlo[5] = xlo[5] - fuzz2;
1624 yn = y + (xlo[5]-x)*v/u;
1625 if (yn>ylo[5]-fuzz && yn<yhi[5]+fuzz) {
1626 orad = hypot(xlo[5]-x,yn-y);
1627 return (orad);
1628 }
1629 if (v>0.0) {
1630 yhi[5] = yhi[5]+fuzz;
1631 xn = x + (yhi[5]-y)*u/v;
1632 if (xn>xlo[5]-fuzz && xn<xhi[5]+fuzz2) {
1633 orad = hypot(xn-x, yhi[5]-y);
1634 return (orad);
1635 }
1636 }
1637 if (v<0.0) {
1638 ylo[5] = ylo[5]-fuzz;
1639 xn = x + (ylo[5]-y)*u/v;
1640 if (xn>xlo[5]-fuzz && xn<xhi[5]+fuzz2) {
1641 orad = hypot(xn-x, ylo[5]-y);
1642 return (orad);
1643 }
1644 }
1645 }
1646 if (xlo[0]>-990.0) {
1647 xlo[0] = xlo[0] - fuzz2;
1648 yn = y + (xlo[0]-x)*v/u;
1649 if (yn>ylo[0]-fuzz && yn<yhi[0]+fuzz) {
1650 orad = hypot(xlo[0]-x,yn-y);
1651 return (orad);
1652 }
1653 if (v>0.0) {
1654 yhi[0] = yhi[0]+fuzz;
1655 xn = x + (yhi[0]-y)*u/v;
1656 if (xn>xlo[0]-fuzz && xn<xhi[0]+fuzz2) {
1657 orad = hypot(xn-x, yhi[0]-y);
1658 return (orad);
1659 }
1660 }
1661 if (v<0.0) {
1662 ylo[0] = ylo[0]-fuzz;
1663 xn = x + (ylo[0]-y)*u/v;
1664 if (xn>xlo[0]-fuzz && xn<xhi[0]+fuzz2) {
1665 orad = hypot(xn-x, ylo[0]-y);
1666 return (orad);
1667 }
1668 }
1669 }
1670 return (orad);
1671 }
1672 if (u>0.0) {
1673 if (xlo[4]>-990.0) {
1674 xhi[4] = xhi[4] + fuzz2;
1675 yn = y + (xhi[4]-x)*v/u;
1676 if (yn>ylo[4]-fuzz && yn<=yhi[4]+fuzz) {
1677 orad = hypot(xhi[4]-x,yn-y);
1678 return (orad);
1679 }
1680 }
1681 if (v>0.0 && xlo[3]>-990.0) {
1682 yhi[3] = yhi[3] + fuzz;
1683 xn = x + (yhi[3]-y)*u/v;
1684 if (xn>xlo[3]-fuzz2 && xn<xhi[3]+fuzz) {
1685 orad = hypot(xn-x, yhi[3]-y);
1686 return (orad);
1687 }
1688 xhi[3] = xhi[3] + fuzz2;
1689 yn = y + (xhi[3]-x)*v/u;
1690 if (yn>ylo[3]-fuzz && yn<yhi[3]+fuzz) {
1691 orad = hypot(xhi[3]-x,yn-y);
1692 return (orad);
1693 }
1694 }
1695 if (v<0.0 && xlo[6]>-990.0) {
1696 ylo[6] = ylo[6] - fuzz;
1697 xn = x + (ylo[6]-y)*u/v;
1698 if (xn>xlo[6]-fuzz2 && xn<xhi[6]+fuzz) {
1699 orad = hypot(xn-x, ylo[6]-y);
1700 return (orad);
1701 }
1702 xhi[6] = xhi[6] + fuzz2;
1703 yn = y + (xhi[6]-x)*v/u;
1704 if (yn>ylo[6]-fuzz && yn<yhi[6]+fuzz) {
1705 orad = hypot(xhi[6]-x,yn-y);
1706 return (orad);
1707 }
1708 }
1709 }
1710 return (orad);
1711
1712 }
1713
1714 /* Plot a time series of one station (for now). Just plot
1715 the first station encountered. */
1716
gatser(struct gacmn * pcm)1717 void gatser (struct gacmn *pcm) {
1718 struct gastn *stn, *stn2;
1719 struct garpt *rpt, *rpt2, *frpt;
1720 int ipen,axmov,im,i,hemflg;
1721 float rmin,rmax,cmin,cmax,cint;
1722 float x1,x2,x,y,tsav,dir;
1723
1724 stn = pcm->result[0].stn;
1725 rpt = stn->rpt;
1726 frpt = rpt; /* Plot stnid = 1st report */
1727 if (rpt == NULL) {
1728 gaprnt (0,"No stations to plot\n");
1729 return;
1730 }
1731
1732 /* First pass just get max and min values for this station */
1733
1734 rmin = 9.99e33;
1735 rmax = -9.99e33;
1736 while (rpt!=NULL) {
1737 if (!cmpch(frpt->stid,rpt->stid,8)) {
1738 if (rpt->val != stn->undef) {
1739 if (rpt->val < rmin) rmin = rpt->val;
1740 if (rpt->val > rmax) rmax = rpt->val;
1741 }
1742 }
1743 rpt = rpt->rpt;
1744 }
1745 if (rmin>rmax) {
1746 printf ("All stations undefined for this variable\n");
1747 return;
1748 }
1749
1750 /* Do scaling */
1751
1752 i = pcm->frame;
1753 if (pcm->tser) pcm->frame = 0;
1754 gas1d (pcm, rmin, rmax, 3, pcm->rotate, NULL, stn);
1755 pcm->frame = i;
1756
1757 /* Set up graphics */
1758
1759 if (pcm->ccolor<0) pcm->ccolor=1;
1760 gxcolr (pcm->ccolor);
1761 if (pcm->cstyle>0) gxstyl(pcm->cstyle);
1762 else gxstyl(1);
1763 gxwide (pcm->cthick);
1764
1765 /* Next pass plot lines. */
1766
1767 if (pcm->tser == 1) {
1768 rpt = stn->rpt;
1769 while (rpt!=NULL) {
1770 gxconv (rpt->tim,0,&x,&y,3);
1771 if (rpt->val != stn->undef) {
1772 i = (int)(rpt->val+0.5);
1773 wxsym (i, x, pcm->ysiz1, pcm->digsiz*1.5, -1, pcm->wxcols);
1774 } else {
1775 gxchpl ("M",1,x,pcm->ysiz1,pcm->digsiz,pcm->digsiz,0.0);
1776 }
1777 rpt = rpt->rpt;
1778 }
1779 } else if (pcm->tser==2) {
1780 stn2 = pcm->result[1].stn;
1781 rpt = stn->rpt;
1782 rpt2 = stn2->rpt;
1783 while (rpt!=NULL) {
1784 if (rpt->val != stn->undef && rpt2->val!=stn2->undef) {
1785 gxconv (rpt->tim,rpt->val,&x,&y,3);
1786 if (rpt->val==0.0 && rpt2->val==0.0) dir = 0.0;
1787 else dir = atan2(rpt2->val,rpt->val);
1788 hemflg = 0;
1789 if (pcm->hemflg == 1) hemflg = 1;
1790 else if (pcm->hemflg == 0) hemflg = 0;
1791 else if (rpt->lat<0.0) hemflg = 1;
1792 gabarb (x, pcm->ysiz1, pcm->digsiz*3.5, pcm->digsiz*2.0,
1793 pcm->digsiz*0.25, dir, hypot(rpt->val,rpt2->val), hemflg);
1794 }
1795 rpt = rpt->rpt;
1796 rpt2 = rpt2->rpt;
1797 }
1798 } else {
1799 gxclip (pcm->xsiz1-0.01, pcm->xsiz2+0.01, pcm->ysiz1, pcm->ysiz2);
1800 if (pcm->cstyle!=0) {
1801 rpt = stn->rpt;
1802 tsav = rpt->tim;
1803 ipen = 3;
1804 while (rpt!=NULL) {
1805 if (!cmpch(frpt->stid,rpt->stid,8)) {
1806 if (rpt->val != stn->undef) {
1807 if (rpt->tim - tsav > 1.0 && !pcm->miconn) ipen = 3;
1808 if (pcm->rotate) gxconv (rpt->val,rpt->tim,&x,&y,3);
1809 else gxconv (rpt->tim,rpt->val,&x,&y,3);
1810 gxplot (x,y,ipen);
1811 ipen = 2;
1812 } else if (!pcm->miconn) ipen=3;
1813 tsav = rpt->tim;
1814 }
1815 rpt = rpt->rpt;
1816 }
1817 }
1818
1819 rpt = stn->rpt;
1820 im = pcm->cmark;
1821 if (im>0) {
1822 while (rpt!=NULL) {
1823 if (!cmpch(frpt->stid,rpt->stid,8)) {
1824 if (rpt->val != stn->undef) {
1825 if (pcm->rotate) gxconv (rpt->val,rpt->tim,&x,&y,3);
1826 else gxconv (rpt->tim,rpt->val,&x,&y,3);
1827 if (im==1 || im==2 || im==4) {
1828 gxcolr (gxqbck());
1829 if (im==1) gxmark (4,x,y,pcm->digsiz+0.01);
1830 else gxmark (im+1,x,y,pcm->digsiz+0.01);
1831 gxcolr(pcm->ccolor);
1832 }
1833 gxmark (im,x,y,pcm->digsiz+0.01);
1834 }
1835 }
1836 rpt = rpt->rpt;
1837 }
1838 }
1839
1840 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
1841 if (pcm->rotate) gaaxpl(pcm,4,3);
1842 else gaaxpl(pcm,3,4);
1843 }
1844 gxwide (pcm->annthk-3);
1845 gxcolr (pcm->anncol);
1846 if (pcm->pass==0 && pcm->grdsflg)
1847 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
1848 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
1849 gagsav (14,pcm,NULL);
1850 }
1851
1852 /* Plot a vertical profile given one or more stations and
1853 a dimension environment where only z varies */
1854
gapprf(struct gacmn * pcm)1855 void gapprf (struct gacmn *pcm) {
1856 struct gastn *stn;
1857 struct garpt *anch, **prev, *next, *rpt, *srpt;
1858 int flag,i,ipen;
1859 float rmin,rmax,cmin,cmax,cint,x,y;
1860 char stid[10];
1861
1862 /* Reorder the linked list of reports so they are sorted */
1863
1864 stn = pcm->result[0].stn;
1865 rpt = stn->rpt;
1866 anch = NULL;
1867 while (rpt!=NULL) {
1868 prev = &anch;
1869 srpt = anch;
1870 flag = 0;
1871 while (srpt!=NULL) {
1872 if (!cmpch(srpt->stid,rpt->stid,8)) {
1873 flag = 1;
1874 if (srpt->lev < rpt->lev) break;
1875 } else if (flag) break;
1876 prev = &(srpt->rpt);
1877 srpt = srpt->rpt;
1878 }
1879 next = rpt->rpt;
1880 rpt->rpt = srpt;
1881 *prev = rpt;
1882 rpt = next;
1883 }
1884 stn->rpt = anch;
1885
1886 if (stn->rpt==NULL) {
1887 gaprnt (0,"No station values to plot\n");
1888 return;
1889 }
1890
1891 /* Get max and min */
1892
1893 rmin = 9.99e33;
1894 rmax = -9.99e33;
1895 rpt = stn->rpt;
1896 while (rpt!=NULL) {
1897 if (rpt->val!=stn->undef) {
1898 if (rmin>rpt->val) rmin = rpt->val;
1899 if (rmax<rpt->val) rmax = rpt->val;
1900 }
1901 rpt=rpt->rpt;
1902 }
1903 if (rmin>rmax) {
1904 gaprnt (0,"All stations values are undefined\n");
1905 return;
1906 }
1907
1908 /* Scaling */
1909
1910 gas1d (pcm, rmin, rmax, 2, !(pcm->rotate), NULL, stn);
1911
1912 /* Do plot */
1913
1914 rpt = stn->rpt;
1915 for (i=0; i<8; i++) stid[i]=' ';
1916 gxstyl(pcm->cstyle);
1917 if (pcm->ccolor<1) gxcolr(1);
1918 else gxcolr(pcm->ccolor);
1919 gxwide (pcm->cthick);
1920 ipen = 3;
1921 while (rpt!=NULL) {
1922 if (rpt->val==stn->undef) {
1923 if (!pcm->miconn) ipen = 3;
1924 } else {
1925 if (pcm->rotate) gxconv (rpt->lev,rpt->val,&x,&y,2);
1926 else gxconv (rpt->val,rpt->lev,&x,&y,2);
1927 if (cmpch(stid,rpt->stid,8)) ipen = 3;
1928 gxplot (x,y,ipen);
1929 ipen=2;
1930 }
1931 for (i=0; i<8; i++) stid[i]=rpt->stid[i];
1932 rpt=rpt->rpt;
1933 }
1934 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
1935 if (pcm->rotate) gaaxpl(pcm,2,4);
1936 else gaaxpl(pcm,4,2);
1937 gxcolr(pcm->anncol);
1938 gxwide (pcm->annthk-3);
1939 if (pcm->pass==0 && pcm->grdsflg)
1940 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
1941 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
1942 gagsav (13,pcm,NULL);
1943 }
1944
1945 /* Routine to set up scaling for a 1-D plot. */
1946
gas1d(struct gacmn * pcm,float cmin,float cmax,int dim,int rotflg,struct gagrid * pgr,struct gastn * stn)1947 void gas1d (struct gacmn *pcm, float cmin, float cmax, int dim,
1948 int rotflg, struct gagrid *pgr, struct gastn *stn) {
1949 float cint,cmn,cmx,x1,x2,xt1,xt2,yt1,yt2;
1950 int axmov;
1951
1952 idiv = 1; jdiv = 1; /* No grid expansion factor */
1953 gxrset(3); /* Reset all scaling */
1954
1955 /* Set plot area limits */
1956
1957 if (pcm->paflg) {
1958 pcm->xsiz1 = pcm->pxmin;
1959 pcm->xsiz2 = pcm->pxmax;
1960 pcm->ysiz1 = pcm->pymin;
1961 pcm->ysiz2 = pcm->pymax;
1962 } else {
1963 if (rotflg) {
1964 pcm->xsiz1 = pcm->xsiz/2.0 - pcm->ysiz/3.0;
1965 pcm->xsiz2 = pcm->xsiz/2.0 + pcm->ysiz/3.0;
1966 if (pcm->xsiz1<1.0) pcm->xsiz1 = 1.0;
1967 if (pcm->xsiz2>pcm->xsiz-0.5) pcm->xsiz2 = pcm->xsiz-0.5;
1968 pcm->ysiz1 = 0.75;
1969 pcm->ysiz2 = pcm->ysiz-0.75;
1970 } else {
1971 pcm->xsiz1 = 2.0;
1972 pcm->xsiz2 = pcm->xsiz - 0.5;
1973 pcm->ysiz1 = 0.75;
1974 pcm->ysiz2 = pcm->ysiz - 0.75;
1975 }
1976 }
1977
1978 /* Determine axis limits and label interval. Use user
1979 supplied AXLIM if available. Also try to use most recent
1980 limits if frame hasn't been cleared -- if rotated, force
1981 use of most recent limits, since we don't plot new axis. */
1982
1983 cint = 0.0;
1984 gacsel (cmin,cmax,&cint,&cmn,&cmx);
1985 if (cint==0.0) {
1986 cmn = cmin-5.0;
1987 cmx = cmin+5.0;
1988 cint = 2.0;
1989 } else {
1990 cmn = cmn - 2.0*cint;
1991 cmx = cmx + 2.0*cint;
1992 }
1993 axmov = 1;
1994 /*
1995 * printf("qqqP %g %g %g %g %d %d %d\n",cmin,cmax,cmn,cmx,rotflg,pcm->pass,pcm->aflag);
1996 *printf("qqqP rmin rmax rint: %g %g %g\n",pcm->rmin,pcm->rmax,pcm->rint);
1997 */
1998
1999 if (pcm->aflag == -1 || (rotflg && pcm->pass>0)) {
2000
2001 /*mf 970416
2002 *
2003 * put check for 0 value if doing a second + plot on the same page
2004 * this means that the scaling from the first plot is not appropriate
2005 *
2006 mf*/
2007 if(!(pcm->rmin == 0 && pcm->rmax == 0)) {
2008 cmn = pcm->rmin;
2009 cmx = pcm->rmax;
2010 cint = pcm->rint;
2011 }
2012 axmov = 0;
2013
2014 } else if (pcm->aflag == 1) {
2015 if ( (cmin > (pcm->rmin-pcm->rint*2.0)) &&
2016 (cmax < (pcm->rmax+pcm->rint*2.0)) ) {
2017 cmn = pcm->rmin;
2018 cmx = pcm->rmax;
2019 cint = pcm->rint;
2020 axmov = 0;
2021 }
2022 }
2023 /*
2024 * printf("qqqP after cmn cmx cint: %g %g %g\n",cmn,cmx,cint);
2025 */
2026
2027 if (!pcm->ylpflg && axmov && pcm->yllow>0.0) {
2028 pcm->ylpos = pcm->ylpos - pcm->yllow - 0.1;
2029 }
2030 pcm->yllow = 0.0;
2031
2032 /* Set absolute coordinate scaling for this grid.
2033 Note that this code assumes knowledge of the time
2034 coordinate scaling setup -- namely, that the same
2035 vals are used for t2gr as gr2t. */
2036
2037 if (pgr!=NULL) {
2038 if (dim==3) {
2039 x1 = t2gr(pgr->ivals,&(pcm->tmin));
2040 x2 = t2gr(pgr->ivals,&(pcm->tmax));
2041 /* x1 = x1 + 1 - (float)(pgr->dimmin[3]);
2042 x2 = x2 + 1 - (float)(pgr->dimmin[3]); */
2043 } else {
2044 x1 = pcm->dmin[dim];
2045 x2 = pcm->dmax[dim];
2046 if (dim==2 && pcm->zlog) {
2047 if (x1<=0.0 || x2<=0.0) {
2048 gaprnt (1,"Cannot use log scaling when coordinates <= 0\n");
2049 gaprnt (1,"Linear scaling used\n");
2050 } else {
2051 if (rotflg) {
2052 gxproj(galny);
2053 gxback(gaalny);
2054 } else {
2055 gxproj(galnx);
2056 gxback(gaalnx);
2057 }
2058 x1 = log(x1);
2059 x2 = log(x2);
2060 }
2061 }
2062 if (dim==1 && pcm->coslat) {
2063 if (x1 < -90.0 || x2 > 90.0) {
2064 gaprnt (1,"Cannot use cos lat scaling when coordinates exceed -90 to 90\n");
2065 gaprnt (1,"Linear scaling used\n");
2066 } else {
2067 if (rotflg) {
2068 gxproj(gacly);
2069 gxback(gaacly);
2070 } else {
2071 gxproj(gaclx);
2072 gxback(gaaclx);
2073 }
2074 x1 = sin(x1*3.1416/180.0);
2075 x2 = sin(x2*3.1416/180.0);
2076 }
2077 }
2078 }
2079 if (rotflg) {
2080 pcm->xdim = 4;
2081 pcm->ydim = dim;
2082 pcm->ygr2ab = pgr->igrab;
2083 pcm->yab2gr = pgr->iabgr;
2084 pcm->ygrval = pgr->ivals;
2085 pcm->yabval = pgr->iavals;
2086 xt1 = cmn; xt2 = cmx; yt1 = x1; yt2 = x2;
2087 if (pcm->xflip) {xt1 = cmx; xt2 = cmn;}
2088 if (pcm->yflip) {yt1 = x2; yt2 = x1;}
2089 } else {
2090 pcm->xdim = dim;
2091 pcm->ydim = 4;
2092 pcm->xgr2ab = pgr->igrab;
2093 pcm->xab2gr = pgr->iabgr;
2094 pcm->xgrval = pgr->ivals;
2095 pcm->xabval = pgr->iavals;
2096 xt1 = x1; xt2 = x2; yt1 = cmn; yt2 = cmx;
2097 if (pcm->xflip) {xt1 = x2; xt2 = x1;}
2098 if (pcm->yflip) {yt1 = cmx; yt2 = cmn;}
2099 }
2100 gxscal (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2,xt1,xt2,yt1,yt2);
2101 if (rotflg) {
2102 if (dim==3) jconv = NULL;
2103 else {
2104 jconv = pgr->igrab;
2105 jvars = pgr->ivals;
2106 }
2107 joffset = pgr->dimmin[dim]-1;
2108 iconv = NULL;
2109 ioffset = 0;
2110 } else {
2111 if (dim==3) iconv = NULL;
2112 else {
2113 iconv = pgr->igrab;
2114 ivars = pgr->ivals;
2115 }
2116 ioffset = pgr->dimmin[dim]-1;
2117 jconv = NULL;
2118 joffset = 0;
2119 }
2120 gxgrid (gaconv);
2121 } else {
2122 if (dim==3) {
2123 x1 = t2gr(stn->tvals, &(pcm->tmin));
2124 x2 = t2gr(stn->tvals, &(pcm->tmax));
2125 } else {
2126 x1 = pcm->dmin[dim];
2127 x2 = pcm->dmax[dim];
2128 if (dim==2 && pcm->zlog) {
2129 if (x1<=0.0 || x2<=0.0) {
2130 gaprnt (1,"Cannot use log scaling when coordinates <= 0\n");
2131 gaprnt (1,"Linear scaling used\n");
2132 } else {
2133 if (rotflg) {
2134 gxproj(galny);
2135 gxback(gaalny);
2136 } else {
2137 gxproj(galnx);
2138 gxback(gaalnx);
2139 }
2140 x1 = log(x1);
2141 x2 = log(x2);
2142 }
2143 }
2144 if (dim==1 && pcm->coslat) {
2145 if (x1 < -90.0 || x2 > 90.0) {
2146 gaprnt (1,"Cannot use cos lat scaling when coordinates exceed -90 to 90\n");
2147 gaprnt (1,"Linear scaling used\n");
2148 } else {
2149 if (rotflg) {
2150 gxproj(gacly);
2151 gxback(gaacly);
2152 } else {
2153 gxproj(gaclx);
2154 gxback(gaaclx);
2155 }
2156 x1 = sin(x1*3.1416/180.0);
2157 x2 = sin(x2*3.1416/180.0);
2158 }
2159 }
2160 }
2161 if (rotflg) {
2162 pcm->xdim = 4;
2163 pcm->ydim = dim;
2164 if (dim==3) pcm->ygrval = stn->tvals;
2165 xt1 = cmn; xt2 = cmx; yt1 = x1; yt2 = x2;
2166 if (pcm->xflip) {xt1 = cmx; xt2 = cmn;}
2167 if (pcm->yflip) {yt1 = x2; yt2 = x1;}
2168 gxscal (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2,xt1,xt2,yt1,yt2);
2169 } else {
2170 pcm->xdim = dim;
2171 pcm->ydim = 4;
2172 if (dim==3) pcm->xgrval = stn->tvals;
2173 xt1 = x1; xt2 = x2; yt1 = cmn; yt2 = cmx;
2174 if (pcm->xflip) {xt1 = x1; xt2 = x2;}
2175 if (pcm->yflip) {yt1 = cmn; yt2 = cmx;}
2176 gxscal (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2,x1,x2,cmn,cmx);
2177 }
2178 }
2179 pcm->rmin = cmn;
2180 pcm->rmax = cmx;
2181 pcm->rint = cint;
2182 if (pcm->aflag == 0) pcm->aflag = 1;
2183 gafram (pcm);
2184 }
2185
2186 /* Set up grid level scaling for a 2-D plot */
2187
gas2d(struct gacmn * pcm,struct gagrid * pgr,int imap)2188 void gas2d (struct gacmn *pcm, struct gagrid *pgr, int imap) {
2189 float x1,x2,y1,y2,xt1,xt2,yt1,yt2;
2190 int idim,jdim;
2191
2192 gxrset (3); /* Reset all scaling */
2193
2194 /* Set up linear level scaling (level 1) and map level scaling
2195 (level 2). If no map drawn, just do linear level scaling. */
2196
2197 idim = pgr->idim;
2198 jdim = pgr->jdim;
2199 pcm->xdim = idim;
2200 pcm->ydim = jdim;
2201 pcm->xgrval = pgr->ivals; pcm->ygrval = pgr->jvals;
2202 pcm->xabval = pgr->iavals; pcm->yabval = pgr->javals;
2203 pcm->xgr2ab = pgr->igrab; pcm->ygr2ab = pgr->jgrab;
2204 pcm->xab2gr = pgr->iabgr; pcm->yab2gr = pgr->jabgr;
2205 if (idim==0 && jdim==1) {
2206 gamscl (pcm); /* Do map level scaling */
2207 if (imap) gawmap(pcm, 1); /* Draw map if requested */
2208 }
2209 else {
2210 if (pcm->paflg) {
2211 pcm->xsiz1 = pcm->pxmin;
2212 pcm->xsiz2 = pcm->pxmax;
2213 pcm->ysiz1 = pcm->pymin;
2214 pcm->ysiz2 = pcm->pymax;
2215 } else {
2216 pcm->xsiz1 = 1.5;
2217 pcm->xsiz2 = pcm->xsiz-0.5;
2218 pcm->ysiz1 = 1.0;
2219 pcm->ysiz2 = pcm->ysiz-0.5;
2220 }
2221 if (idim==3) {
2222 x1 = t2gr(pgr->ivals,&(pcm->tmin));
2223 x2 = t2gr(pgr->ivals,&(pcm->tmax));
2224 /* x1 = x1 + 1 - (float)(pgr->dimmin[3]);
2225 x2 = x2 + 1 - (float)(pgr->dimmin[3]); */
2226 } else if (idim==4) {
2227 x1 = 1; x2 = pgr->isiz; /* COLL */
2228 } else {
2229 x1 = pcm->dmin[idim];
2230 x2 = pcm->dmax[idim];
2231 if (idim==2 && pcm->zlog) {
2232 if (x1<=0.0 || x2<=0.0) {
2233 gaprnt (1,"Cannot use log scaling when coordinates <= 0\n");
2234 gaprnt (1,"Linear scaling used\n");
2235 } else {
2236 gxproj(galnx);
2237 gxback(gaalnx);
2238 x1 = log(x1);
2239 x2 = log(x2);
2240 }
2241 }
2242 else if (idim==1 && pcm->coslat) { /* can't have zlog and coslat both */
2243 if (x1 < -90.0 || x2 > 90.0) {
2244 gaprnt (1,"Cannot use cos lat scaling when coordinates exceed -90 to 90\n");
2245 gaprnt (1,"Linear scaling used\n");
2246 } else {
2247 gxproj(gaclx);
2248 gxback(gaaclx);
2249 x1 = sin(x1*3.1416/180.0);
2250 x2 = sin(x2*3.1416/180.0);
2251 }
2252 }
2253 }
2254 if (jdim==3) {
2255 y1 = t2gr(pgr->jvals,&(pcm->tmin));
2256 y2 = t2gr(pgr->jvals,&(pcm->tmax));
2257 /* y1 = y1 + 1 - (float)(pgr->dimmin[3]);
2258 y2 = y2 + 1 - (float)(pgr->dimmin[3]); */
2259 } else {
2260 y1 = pcm->dmin[jdim];
2261 y2 = pcm->dmax[jdim];
2262 if (jdim==2 && pcm->zlog) {
2263 if (y1<=0.0 || y2<=0.0) {
2264 gaprnt (1,"Cannot use log scaling when coordinates <= 0\n");
2265 gaprnt (1,"Linear scaling used\n");
2266 } else {
2267 gxproj(galny);
2268 gxback(gaalny);
2269 y1 = log(y1);
2270 y2 = log(y2);
2271 }
2272 }
2273 else if (jdim==1 && pcm->coslat) { /* can't have zlog and coslat both */
2274 if (y1 < -90.0 || y2 > 90.0) {
2275 gaprnt (1,"Cannot use cos lat scaling when coordinates exceed -90 to 90\n");
2276 gaprnt (1,"Linear scaling used\n");
2277 } else {
2278 gxproj(gacly);
2279 gxback(gaacly);
2280 y1 = sin(y1*3.1416/180.0);
2281 y2 = sin(y2*3.1416/180.0);
2282 }
2283 }
2284 }
2285 xt1 = x1; xt2 = x2; yt1 = y1; yt2 = y2;
2286 if (pcm->xflip) { xt1 = x2; xt2 = x1; }
2287 if (pcm->yflip) { yt1 = y2; yt2 = y1; }
2288 gxscal (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2,xt1,xt2,yt1,yt2);
2289 }
2290
2291 /* Now set up level 2 grid scaling done through gaconv */
2292
2293 if (idim==4) ioffset = 0; /* COLL */
2294 else ioffset = pgr->dimmin[idim]-1;
2295 if (idim==3) iconv = NULL;
2296 else {
2297 iconv = pgr->igrab;
2298 ivars = pgr->ivals;
2299 }
2300 if (jdim==4) joffset = 1; /* COLL */
2301 else joffset = pgr->dimmin[jdim]-1;
2302 if (jdim==3) jconv = NULL;
2303 else {
2304 jconv = pgr->jgrab;
2305 jvars = pgr->jvals;
2306 }
2307 gxgrid (gaconv);
2308
2309 if (idim==4) { /* COLL -- predefine axis */
2310 pcm->rmin = 1;
2311 pcm->rmax = (float)(pgr->isiz);
2312 pcm->axmin = 1.0;
2313 pcm->axmax = pcm->rmax;
2314 pcm->axflg = 1;
2315 pcm->axint = 1.0;
2316 }
2317 }
2318
2319 /* Line plot fill. Fill above and below a line plot */
2320
galfil(struct gacmn * pcm)2321 void galfil (struct gacmn *pcm) {
2322 struct gagrid *pgr1, *pgr2;
2323 float *gr1,*gr2,rmin,rmax,*v1,*v2,*u,*xy,uu,vv;
2324 int rotflg,flag,sflag,cnt,abflg,i,ii,colr;
2325
2326 /* Check args */
2327
2328 if (pcm->numgrd<2) {
2329 gaprnt (0,"Error plotting linefill: Only one grid provided\n");
2330 return;
2331 }
2332 pgr1 = pcm->result[0].pgr;
2333 pgr2 = pcm->result[1].pgr;
2334 if (pgr1->idim!=pgr2->idim || pgr1->jdim!=-1 ||
2335 pgr2->jdim!=-1 || gagchk(pgr1,pgr2,pgr1->idim)) {
2336 gaprnt (0,"Error plotting linefill: Invalid grids --");
2337 gaprnt (0," different scaling\n");
2338 return;
2339 }
2340
2341 gamnmx(pgr1);
2342 gamnmx(pgr2);
2343 if (pgr1->rmin==pgr1->undef || pgr2->rmin==pgr2->undef) {
2344 gaprnt (1,"Cannot plot data - all undefined values \n");
2345 return;
2346 }
2347 rmin = pgr1->rmin;
2348 if (rmin>pgr2->rmin) rmin = pgr2->rmin;
2349 rmax = pgr1->rmax;
2350 if (rmax<pgr2->rmax) rmax = pgr2->rmax;
2351
2352 /* Do scaling */
2353
2354 rotflg = 0;
2355 if (pgr1->idim==2) rotflg = 1;
2356 if (pcm->rotate) rotflg = !rotflg;
2357 gas1d (pcm, rmin, rmax, pgr1->idim, rotflg, pgr1, NULL);
2358
2359 gxclip (pcm->xsiz1-0.01, pcm->xsiz2+0.01, pcm->ysiz1, pcm->ysiz2);
2360
2361 /* Allocate some buffer areas */
2362
2363 u = (float *)malloc(sizeof(float)*pgr1->isiz);
2364 v1 = (float *)malloc(sizeof(float)*pgr1->isiz);
2365 v2 = (float *)malloc(sizeof(float)*pgr1->isiz);
2366 xy = (float *)malloc(sizeof(float)*(pgr1->isiz*4+8));
2367 if (u==NULL || v1==NULL || v2==NULL || xy==NULL) {
2368 gaprnt(0,"Memory allocation error in linefill\n");
2369 return;
2370 }
2371
2372 /* Find a start point. It is the first point where
2373 two points in a row are not missing (in both lines)
2374 and are not equal. */
2375
2376 gr1 = pgr1->grid;
2377 gr2 = pgr2->grid;
2378
2379 i = 1;
2380 while (1) {
2381 if (i>=pgr1->isiz) break;
2382 if (*gr1!=pgr1->undef && *gr2!=pgr2->undef &&
2383 *(gr1+1)!=pgr1->undef && *(gr2+1)!=pgr2->undef &&
2384 (*gr1!=*gr2 || *(gr1+1)!=*(gr2+1))) break;
2385 i++; gr1++; gr2++;
2386 }
2387 if (i>=pgr1->isiz) {
2388 free (u);
2389 free (v1);
2390 free (v2);
2391 free (xy);
2392 gaprnt (1,"Cannot plot data - too many undefined values \n");
2393 return;
2394 }
2395
2396 /* Set up for loop */
2397
2398 cnt = 1;
2399 *u = i;
2400 *v1 = *gr1;
2401 *v2 = *gr2;
2402 abflg = 0;
2403 if (*gr1>*gr2) abflg = 1;
2404 if (*gr1<*gr2) abflg = 2;
2405 i++; gr1++; gr2++;
2406 if (abflg==0) { /* if 1st point is same for both lines */
2407 if (*gr1>*gr2) abflg = 1;
2408 if (*gr1<*gr2) abflg = 2;
2409 }
2410
2411 /* Go through rest of data */
2412
2413 while (i<=pgr1->isiz) {
2414 sflag = 0;
2415
2416 /* If we hit missing data, terminate the current
2417 poly fill */
2418
2419 if (*gr1==pgr1->undef || *gr2==pgr2->undef) {
2420 if (abflg==1) colr=pcm->lfc1;
2421 if (abflg==2) colr=pcm->lfc2;
2422 if (colr>-1) {
2423 gxcolr(colr);
2424 lfout (u, v1, v2, xy, cnt, rotflg);
2425 }
2426 sflag = 1;
2427
2428 /* No missing data. Polygon continues until
2429 the curves cross. */
2430
2431 } else {
2432 if (abflg==0) { /* this shouldn't really happen? */
2433 if (*gr1>*gr2) abflg = 1;
2434 else if (*gr1<*gr2) abflg = 2;
2435 else sflag = 1;
2436 printf ("logic error?\n");
2437 } else if (abflg==1) {
2438 if (*gr1<*gr2) { /* Curves crossed */
2439 lfint (*(v1+cnt-1),*gr1,*(v2+cnt-1),*gr2,&uu,&vv);
2440 uu = uu + (float)(i-1);
2441 *(u+cnt) = uu;
2442 *(v1+cnt) = vv;
2443 *(v2+cnt) = vv;
2444 cnt++;
2445 if (abflg==1) colr=pcm->lfc1;
2446 if (abflg==2) colr=pcm->lfc2;
2447 if (colr>-1) {
2448 gxcolr(colr);
2449 lfout (u, v1, v2, xy, cnt, rotflg);
2450 }
2451 *u = uu; *v1 = vv; *v2 = vv;
2452 *(u+1) = i; *(v1+1) = *gr1; *(v2+1) = *gr2;
2453 cnt = 2;
2454 abflg = 2;
2455 } else if (*gr1==*gr2) { /* Curves touching */
2456 *(v1+cnt) = *gr1;
2457 *(v2+cnt) = *gr2;
2458 *(u+cnt) = i;
2459 cnt++;
2460 if (abflg==1) colr=pcm->lfc1;
2461 if (abflg==2) colr=pcm->lfc2;
2462 if (colr>-1) {
2463 gxcolr(colr);
2464 lfout (u, v1, v2, xy, cnt, rotflg);
2465 }
2466 sflag = 1;
2467 } else { /* curves not crossing; add to poygon */
2468 *(u+cnt) = i;
2469 *(v1+cnt) = *gr1;
2470 *(v2+cnt) = *gr2;
2471 cnt++;
2472 }
2473 } else if (abflg==2) {
2474 if (*gr1>*gr2) {
2475 lfint (*(v1+cnt-1),*gr1,*(v2+cnt-1),*gr2,&uu,&vv);
2476 uu = uu + (float)(i-1);
2477 *(u+cnt) = uu;
2478 *(v1+cnt) = vv;
2479 *(v2+cnt) = vv;
2480 cnt++;
2481 if (abflg==1) colr=pcm->lfc1;
2482 if (abflg==2) colr=pcm->lfc2;
2483 if (colr>-1) {
2484 gxcolr(colr);
2485 lfout (u, v1, v2, xy, cnt, rotflg);
2486 }
2487 *u = uu; *v1 = vv; *v2 = vv;
2488 *(u+1) = i; *(v1+1) = *gr1; *(v2+1) = *gr2;
2489 cnt = 2;
2490 abflg = 1;
2491 } else if (*gr1==*gr2) {
2492 *(v1+cnt) = *gr1;
2493 *(v2+cnt) = *gr2;
2494 *(u+cnt) = i;
2495 cnt++;
2496 if (abflg==1) colr=pcm->lfc1;
2497 if (abflg==2) colr=pcm->lfc2;
2498 if (colr>-1) {
2499 gxcolr(colr);
2500 lfout (u, v1, v2, xy, cnt, rotflg);
2501 }
2502 sflag = 1;
2503 } else {
2504 *(u+cnt) = i;
2505 *(v1+cnt) = *gr1;
2506 *(v2+cnt) = *gr2;
2507 cnt++;
2508 }
2509 }
2510 }
2511
2512 /* If necessary, search for new start point */
2513
2514 if (sflag) {
2515 while (1) {
2516 if (i>=pgr1->isiz) break;
2517 if (*gr1!=pgr1->undef && *gr2!=pgr2->undef &&
2518 *(gr1+1)!=pgr1->undef && *(gr2+1)!=pgr2->undef &&
2519 (*gr1!=*gr2 || *(gr1+1)!=*(gr2+1))) break;
2520 i++; gr1++; gr2++;
2521 }
2522 if (i<pgr1->isiz) {
2523 cnt = 1;
2524 *u = i;
2525 *v1 = *gr1;
2526 *v2 = *gr2;
2527 abflg = 0;
2528 if (*gr1>*gr2) abflg = 1;
2529 if (*gr1<*gr2) abflg = 2;
2530 i++; gr1++; gr2++;
2531 if (abflg==0) {
2532 if (*gr1>*gr2) abflg = 1;
2533 if (*gr1<*gr2) abflg = 2;
2534 }
2535 } else {
2536 cnt = 0;
2537 i = pgr1->isiz+1;
2538 }
2539 } else {
2540 i++; gr1++; gr2++;
2541 }
2542 }
2543
2544 /* Finish up and plot last poly */
2545
2546 if (cnt>1) {
2547 if (abflg==1) colr=pcm->lfc1;
2548 if (abflg==2) colr=pcm->lfc2;
2549 if (colr>-1) {
2550 gxcolr(colr);
2551 lfout (u, v1, v2, xy, cnt, rotflg);
2552 }
2553 }
2554 free (u);
2555 free (v1);
2556 free (v2);
2557 free (xy);
2558 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
2559 if (rotflg) gaaxpl(pcm,4,pgr1->idim);
2560 else gaaxpl(pcm,pgr1->idim,4);
2561 gxwide (pcm->annthk-3);
2562 gxcolr (pcm->anncol);
2563 if (pcm->pass==0 && pcm->grdsflg)
2564 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
2565 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
2566 gagsav (17,pcm,pgr1);
2567 }
2568
lfint(float v11,float v12,float v21,float v22,float * u,float * v)2569 void lfint (float v11, float v12, float v21, float v22,
2570 float *u, float *v) {
2571 *u = (v21-v11)/(v12-v11-v22+v21);
2572 *v = v11 + *u*(v12-v11);
2573 }
2574
lfout(float * u,float * v1,float * v2,float * xyb,int cnt,int rotflg)2575 void lfout (float *u, float *v1, float *v2, float *xyb, int cnt, int rotflg) {
2576 int i,j,ii;
2577 float *xy;
2578
2579 xy = xyb;
2580 if (cnt<2) {
2581 printf ("Internal Logic check 4 in linefill\n");
2582 return;
2583 }
2584 j = 0;
2585 for (i=0; i<cnt; i++) {
2586 if (rotflg) gxconv (*(v1+i),*(u+i),xy,xy+1,3);
2587 else gxconv (*(u+i),*(v1+i),xy,xy+1,3);
2588 j++; xy+=2;
2589 }
2590 i = cnt-1;
2591 if (*(v1+i) != *(v2+i) ) {
2592 if (rotflg) gxconv (*(v2+i),*(u+i),xy,xy+1,3);
2593 else gxconv (*(u+i),*(v2+i),xy,xy+1,3);
2594 j++; xy+=2;
2595 }
2596 for (i=2; i<cnt; i++) {
2597 ii = cnt-i;
2598 if (rotflg) gxconv (*(v2+ii),*(u+ii),xy,xy+1,3);
2599 else gxconv (*(u+ii),*(v2+ii),xy,xy+1,3);
2600 j++; xy+=2;
2601 }
2602 if (*v1 != *v2) {
2603 if (rotflg) gxconv (*v2,*u,xy,xy+1,3);
2604 else gxconv (*u,*v2,xy,xy+1,3);
2605 j++; xy+=2;
2606 }
2607 gxfill (xyb,j);
2608 }
2609
2610 /* Do line graph (barflg=0) or bar graph (barflg=1) or
2611 1-D wind vector/barb plot (barflg=2) */
2612
gagrph(struct gacmn * pcm,int barflg)2613 void gagrph (struct gacmn *pcm, int barflg) {
2614 struct gagrid *pgr,*pgr2,*pgr3;
2615 float cmin,cmax,cint,x1,x2,y1,y2,x,y,xz,yz,rmin,rmax;
2616 float *gr,*gr2,*gr3,gap,siz,umax,vmax,vscal,dir,xx,yy;
2617 int idim,ip,i,im,axmov,bflg,rotflg,flag,pflg,hflg,hemflg;
2618
2619 /* Determine min and max */
2620
2621 pgr = pcm->result[0].pgr;
2622 gamnmx(pgr);
2623 if (pgr->rmin==pgr->undef) {
2624 gaprnt (1,"Cannot plot data - all undefined values \n");
2625 gxcolr(1);
2626 if (pcm->dwrnflg) gxchpl ("Entire Grid Undefined",21,3.0,4.5,0.3,0.25,0.0);
2627 return;
2628 }
2629 rmin = pgr->rmin; rmax = pgr->rmax;
2630
2631 flag = 0;
2632 if (pcm->numgrd>1) {
2633 flag = 1;
2634 pgr2 = pcm->result[1].pgr;
2635 if (pgr2->idim!=pgr->idim || pgr2->jdim!=pgr->jdim ||
2636 gagchk(pgr2,pgr,pgr->idim) || gagchk(pgr2,pgr,pgr->jdim)) {
2637 flag = 0;
2638 gaprnt (0,"Error in line/bar graph: Invalid 2nd field");
2639 gaprnt (0," 2nd grid ignored -- has different scaling");
2640 } else {
2641 if (barflg==1) {
2642 gamnmx(pgr2);
2643 if (pgr2->rmin==pgr2->undef) {
2644 gaprnt (1,"Cannot plot data - 2nd arg all undefined values \n");
2645 return;
2646 }
2647 if (rmin>pgr2->rmin) rmin = pgr2->rmin;
2648 if (rmax<pgr2->rmax) rmax = pgr2->rmax;
2649 }
2650 }
2651 if (pcm->numgrd>2 && flag==1) {
2652 flag = 2;
2653 pgr3 = pcm->result[2].pgr;
2654 if (pgr3->idim!=pgr->idim || pgr3->jdim!=pgr->jdim ||
2655 gagchk(pgr3,pgr,pgr->idim) || gagchk(pgr3,pgr,pgr->jdim)) {
2656 flag = 1;
2657 gaprnt (0,"Error in line/bar graph: Invalid 3rd field");
2658 gaprnt (0," 3rd grid ignored -- has different scaling");
2659 }
2660 }
2661 }
2662
2663 /* Do scaling */
2664
2665 rotflg = 0;
2666 if (pgr->idim==2) rotflg = 1;
2667 if (pcm->rotate) rotflg = !rotflg;
2668
2669 gas1d (pcm, rmin, rmax, pgr->idim, rotflg, pgr, NULL);
2670
2671 /* Set up graphics */
2672
2673 if (pcm->ccolor<0) pcm->ccolor=1;
2674 gxcolr (pcm->ccolor);
2675 gxstyl(pcm->cstyle);
2676 gxwide (pcm->cthick);
2677 gr=pgr->grid;
2678 if (flag) gr2 = pgr2->grid;
2679 gxclip (pcm->xsiz1-0.01, pcm->xsiz2+0.01, pcm->ysiz1, pcm->ysiz2);
2680
2681 /* Do bar graph */
2682
2683 if (barflg) {
2684 bflg = pcm->barflg;
2685 gap = ((float)pcm->bargap)*0.5/100.0;
2686 gap = 0.5-gap;
2687 if (rotflg) {
2688 if (bflg==1) gxconv (pcm->barbase,1.0,&xz,&y,3);
2689 else if (bflg==0) xz = pcm->xsiz1;
2690 else xz = pcm->xsiz2;
2691 if (bflg==1 && (xz<pcm->xsiz1||xz>pcm->xsiz2)) {
2692 gaprnt(0,"Error drawing bargraph: base value out of range\n");
2693 gaprnt(0," Drawing graph from bottom\n");
2694 printf ("%g \n",xz);
2695 bflg = 0;
2696 xz = pcm->xsiz1;
2697 }
2698 for (i=1;i<=pgr->isiz;i++) {
2699 pflg = 1;
2700 if (flag) {
2701 if (*gr == pgr->undef || *gr2 == pgr2->undef) {
2702 pflg=0;
2703 } else {
2704 gxconv (*gr,(float)(i)-gap,&x2,&y1,3);
2705 gxconv (*gr,(float)(i)+gap,&x2,&y2,3);
2706 gxconv (*gr,(float)(i),&x2,&yz,3);
2707 gxconv (*gr2,(float)(i),&x1,&yz,3);
2708 }
2709 } else {
2710 if (*gr == pgr->undef) {
2711 pflg = 0;
2712 } else {
2713 gxconv (*gr,(float)(i)-gap,&x2,&y1,3);
2714 gxconv (*gr,(float)(i)+gap,&x2,&y2,3);
2715 gxconv (*gr,(float)(i),&x2,&yz,3);
2716 x1 = xz;
2717 }
2718 }
2719 if (pflg) {
2720 if (barflg==2) {
2721 gxplot(x1,yz,3);
2722 gxplot(x2,yz,2);
2723 gxplot(x1,y1,3);
2724 gxplot(x1,y2,2);
2725 gxplot(x2,y1,3);
2726 gxplot(x2,y2,2);
2727 } else if (pcm->barolin) {
2728 gxplot(x1,y1,3);
2729 gxplot(x1,y2,2);
2730 gxplot(x2,y2,2);
2731 gxplot(x2,y1,2);
2732 gxplot(x1,y1,2);
2733 } else gxrecf (xz, x, y1, y2);
2734 }
2735 /*
2736 } else {
2737 gxconv (0.0,(float)(i)-gap,&x,&y1,3);
2738 gxconv (0.0,(float)(i),&x,&y,3);
2739 gxconv (0.0,(float)(i)+gap,&x,&y2,3);
2740 siz = fabs(y2-y1)*0.7;
2741 gxchpl ("M",1,xz+siz*0.3,y-siz*0.5,siz,siz,0.0);
2742 }
2743 */
2744 gr++;
2745 if (flag) gr2++;
2746 }
2747 } else {
2748 if (bflg==1) gxconv (1.0,pcm->barbase,&x,&yz,3);
2749 else if (bflg==0) yz = pcm->ysiz1;
2750 else yz = pcm->ysiz2;
2751 if (bflg==1 && (yz<pcm->ysiz1||yz>pcm->ysiz2)) {
2752 gaprnt(0,"Error drawing bargraph: base value out of range\n");
2753 gaprnt(0," Drawing graph from bottom\n");
2754 printf ("%g \n",yz);
2755 bflg = 0;
2756 yz = pcm->ysiz1;
2757 }
2758 for (i=1;i<=pgr->isiz;i++) {
2759 pflg = 1;
2760 if (flag) {
2761 if (*gr == pgr->undef || *gr2 == pgr2->undef) {
2762 pflg=0;
2763 } else {
2764 gxconv ((float)(i)-gap,*gr,&x1,&y1,3);
2765 gxconv ((float)(i)+gap,*gr,&x2,&y1,3);
2766 gxconv ((float)(i),*gr,&xz,&y1,3);
2767 gxconv ((float)(i),*gr2,&xz,&y2,3);
2768 }
2769 } else {
2770 if (*gr == pgr->undef) {
2771 pflg = 0;
2772 } else {
2773 gxconv ((float)(i)-gap,*gr,&x1,&y1,3);
2774 gxconv ((float)(i)+gap,*gr,&x2,&y1,3);
2775 gxconv ((float)(i),*gr,&xz,&y1,3);
2776 y2 = yz;
2777 }
2778 }
2779 if (pflg) {
2780 if (barflg==2) {
2781 gxplot(xz,y1,3);
2782 gxplot(xz,y2,2);
2783 gxplot(x1,y1,3);
2784 gxplot(x2,y1,2);
2785 gxplot(x1,y2,3);
2786 gxplot(x2,y2,2);
2787 } else if (pcm->barolin) {
2788 gxplot(x1,y1,3);
2789 gxplot(x1,y2,2);
2790 gxplot(x2,y2,2);
2791 gxplot(x2,y1,2);
2792 gxplot(x1,y1,2);
2793 } else gxrecf (x1, x2, y1, y2);
2794 }
2795 /*
2796 } else {
2797 gxconv ((float)(i)-gap,0.0,&x1,&y,3);
2798 gxconv ((float)(i),0.0,&x,&y,3);
2799 gxconv ((float)(i)+gap,0.0,&x2,&y,3);
2800 siz = fabs(x2-x1)*0.7;
2801 gxchpl ("M",1,x-siz*0.5,yz+siz*0.3,siz,siz,0.0);
2802 }
2803 */
2804 gr++;
2805 if (flag) gr2++;
2806 }
2807 }
2808
2809 /* Do line graph, or wind vectors/barbs when 3 grids */
2810
2811 } else {
2812
2813 /* Draw the line first */
2814
2815 if (pcm->cstyle!=0 && flag<2) {
2816 ip=3;
2817 for (i=1;i<=pgr->isiz;i++) {
2818 if (*gr == pgr->undef) {
2819 if (!pcm->miconn) ip = 3;
2820 } else {
2821 if (rotflg) gxconv (*gr,(float)i,&x,&y,3);
2822 else gxconv ((float)i,*gr,&x,&y,3);
2823 gxplot (x,y,ip);
2824 ip=2;
2825 }
2826 gr++;
2827 }
2828 }
2829
2830 /* Now draw the markers, or wind vectors/barbs */
2831
2832 im = pcm->cmark;
2833 if (im>0 || flag==2) {
2834 gxstyl (1);
2835 gr=pgr->grid;
2836 if (flag==2) { /* if vectors/barbs, get ready */
2837 gr2 = pgr2->grid;
2838 gr3 = pgr3->grid;
2839
2840 /* hflg=0 idim is lat
2841 hflg=2 plot nhem
2842 hflg=3 plot shem */
2843
2844 if (pcm->hemflg==0) hflg = 2;
2845 else if (pcm->hemflg==1) hflg = 3;
2846 else {
2847 if (pgr2->idim==1) hflg = 0;
2848 else {
2849 if (pcm->dmin[1] < 0.0) hflg = 3;
2850 else hflg = 2;
2851 }
2852 }
2853 if (!pcm->arrflg) {
2854 gamnmx (pgr2);
2855 gamnmx (pgr3);
2856 umax = pgr2->rmax;
2857 if (umax<fabs(pgr2->rmin)) umax = fabs(pgr2->rmin);
2858 vmax = pgr3->rmax;
2859 if (vmax<fabs(pgr3->rmin)) vmax = fabs(pgr3->rmin);
2860 vscal = hypot(umax,vmax);
2861 if (vscal>0.0) {
2862 x = floor(log10(vscal));
2863 y = floor(vscal/pow(10.0,x));
2864 vscal = y * pow(10.0,x);
2865 pcm->arrsiz = 0.5;
2866 } else vscal=1.0;
2867 pcm->arrmag = vscal;
2868 } else {
2869 vscal = pcm->arrmag;
2870 }
2871 pcm->arrflg = 1;
2872 }
2873 for (i=1;i<=pgr->isiz;i++) {
2874 if (*gr != pgr->undef) {
2875 if (rotflg) gxconv (*gr,(float)i,&x,&y,3);
2876 else gxconv ((float)i,*gr,&x,&y,3);
2877 if (flag==2) { /* handle vectors/barbs */
2878 /* xcv */
2879 if (*gr2!=pgr2->undef && *gr3!=pgr3->undef) {
2880 if (rotflg) gxgrmp (*gr,(float)i,&xx,&yy);
2881 else gxgrmp ((float)i,*gr,&yy,&xx);
2882 hemflg = 0;
2883 if (hflg==0 && yy<0.0) hemflg = 1;
2884 if (hflg==3) hemflg = 1;
2885 if (*gr2==0.0 && *gr3==0.0) dir = 0.0;
2886 else dir = atan2(*gr3,*gr2);
2887 if (pcm->gout1a==2) {
2888 gabarb (x, y, pcm->digsiz*3.5, pcm->digsiz*2.0,
2889 pcm->digsiz*0.25, dir, hypot(*gr2,*gr3), hemflg);
2890 } else {
2891 if (vscal>0.0) {
2892 gaarrw (x, y, dir, pcm->arrsiz*hypot(*gr2,*gr3)/vscal,
2893 pcm->ahdsiz);
2894 } else {
2895 gaarrw (x, y, dir, pcm->arrsiz, pcm->ahdsiz);
2896 }
2897 }
2898 }
2899 } else {
2900 if (im==1 || im==2 || im==4 || im==8) {
2901 gxcolr (gxqbck());
2902 if (im==1) gxmark (4,x,y,pcm->digsiz+0.01);
2903 else gxmark (im+1,x,y,pcm->digsiz+0.01);
2904 gxcolr(pcm->ccolor);
2905 }
2906 gxmark (im,x,y,pcm->digsiz+0.01);
2907 }
2908 }
2909 gr++;
2910 if (flag==2) { gr2++; gr3++; }
2911 }
2912 }
2913 }
2914 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
2915 if (rotflg) gaaxpl(pcm,4,pgr->idim);
2916 else gaaxpl(pcm,pgr->idim,4);
2917 gxwide (pcm->annthk-3);
2918 gxcolr (pcm->anncol);
2919 if (pcm->pass==0 && pcm->grdsflg)
2920 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
2921 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
2922 if (barflg) gagsav (8,pcm,pgr);
2923 else gagsav (6,pcm,pgr);
2924 }
2925
2926 /* Do 2-D streamlines */
2927
gastrm(struct gacmn * pcm)2928 void gastrm (struct gacmn *pcm) {
2929 struct gagrid *pgru, *pgrv, *pgrc;
2930 int flag,lcol,*cpnt,irb;
2931 float *u, *v, *c;
2932
2933 if (pcm->numgrd<2) {
2934 gaprnt (0,"Error plotting streamlines: Only one grid provided\n");
2935 return;
2936 }
2937 pgru = pcm->result[0].pgr;
2938 pgrv = pcm->result[1].pgr;
2939 if (pgru->idim!=pgrv->idim || pgru->jdim!=pgrv->jdim ||
2940 gagchk(pgru,pgrv,pgru->idim) || gagchk(pgru,pgrv,pgru->jdim)) {
2941 gaprnt (0,"Error plotting streamlines: Invalid grids\n");
2942 gaprnt (0," Vector component grids have difference scaling\n");
2943 return;
2944 }
2945 flag = 0;
2946 if (pcm->numgrd>2) {
2947 flag = 1;
2948 pgrc = pcm->result[2].pgr;
2949 if (pgrc->idim!=pgru->idim || pgrc->jdim!=pgru->jdim ||
2950 gagchk(pgrc,pgru,pgru->idim) || gagchk(pgrc,pgru,pgru->jdim)) {
2951 flag = 0;
2952 gaprnt (0,"Error plotting streamlines: Invalid color grid");
2953 gaprnt (0," Color grid ignored -- has different scaling");
2954 }
2955 }
2956
2957 if ( (pcm->rotate && (pgru->idim!=2 || pgru->jdim!=3)) ||
2958 (!pcm->rotate && pgru->idim==2 && pgru->jdim==3)) {
2959 pgru = gaflip(pgru,pcm);
2960 pgrv = gaflip(pgrv,pcm);
2961 if (flag) pgrc = gaflip(pgrc,pcm);
2962 }
2963
2964 gxstyl (1);
2965 gxwide (1);
2966 gas2d (pcm, pgru, 1); /* Set up scaling, draw map if apprprt */
2967
2968 gafram (pcm);
2969 idiv = 1.0; jdiv = 1.0;
2970
2971 if (flag) {
2972 gamnmx (pgrc);
2973 if (pgrc->rmin==pgrc->undef) {
2974 gaprnt (0,"Connot color vectors -- Color grid all undefined\n");
2975 flag = 0;
2976 } else gaselc(pcm,pgrc->rmin,pgrc->rmax);
2977 }
2978
2979 u = pgru->grid;
2980 v = pgrv->grid;
2981 if (flag) c = pgrc->grid;
2982
2983 if (pcm->ccolor>=0) lcol = pcm->ccolor;
2984 else lcol=1;
2985 gxcolr (lcol);
2986 gxstyl(pcm->cstyle);
2987 gxwide(pcm->cthick);
2988 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
2989
2990 if (flag) {
2991 gxstrm (u,v,c,pgru->isiz,pgru->jsiz,pgru->undef,pgrv->undef,
2992 pgrc->undef,flag,pcm->shdlvs,pcm->shdcls,pcm->shdcnt,pcm->strmden);
2993 } else {
2994 gxstrm (u,v,NULL,pgru->isiz,pgru->jsiz,pgru->undef,pgrv->undef,
2995 0.0,flag,pcm->shdlvs,pcm->shdcls,pcm->shdcnt,pcm->strmden);
2996 }
2997
2998 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
2999 gxwide (pcm->annthk-3);
3000 gxcolr (pcm->anncol);
3001 if (pcm->pass==0 && pcm->grdsflg)
3002 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
3003 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
3004 gaaxpl(pcm,pgru->idim,pgru->jdim);
3005 gagsav (9,pcm,pgru);
3006 }
3007
3008 /* Do 2-D vector plot */
3009
gavect(struct gacmn * pcm,int brbflg)3010 void gavect (struct gacmn *pcm, int brbflg) {
3011 struct gagrid *pgru, *pgrv, *pgrc;
3012 float x, y, vscal, *u, *v, *c, lon, lat, adj, rrb, dir, umax, vmax;
3013 int i, j, flag, lcol, len, ii, irb, hemflg, hflg;
3014
3015 if (pcm->numgrd<2) {
3016 gaprnt (0,"Error plotting vector field: Only one grid provided\n");
3017 return;
3018 }
3019 pgru = pcm->result[0].pgr;
3020 pgrv = pcm->result[1].pgr;
3021 if (pgru->idim!=pgrv->idim || pgru->jdim!=pgrv->jdim ||
3022 gagchk(pgru,pgrv,pgru->idim) || gagchk(pgru,pgrv,pgru->jdim)) {
3023 gaprnt (0,"Error plotting vector/barb field: Invalid grids\n");
3024 gaprnt (0," Vector component grids have difference scaling\n");
3025 return;
3026 }
3027 flag = 0;
3028 if (pcm->numgrd>2) {
3029 flag = 1;
3030 pgrc = pcm->result[2].pgr;
3031 if (pgrc->idim!=pgru->idim || pgrc->jdim!=pgru->jdim ||
3032 gagchk(pgrc,pgru,pgru->idim) || gagchk(pgrc,pgru,pgru->jdim)) {
3033 flag = 0;
3034 gaprnt (0,"Error plotting vector/barb field: Invalid color grid");
3035 gaprnt (0," Color grid ignored -- has different scaling");
3036 }
3037 }
3038
3039 if ( (pcm->rotate && (pgru->idim!=2 || pgru->jdim!=3)) ||
3040 (!pcm->rotate && pgru->idim==2 && pgru->jdim==3)) {
3041 pgru = gaflip(pgru,pcm);
3042 pgrv = gaflip(pgrv,pcm);
3043 if (flag) pgrc = gaflip(pgrc,pcm);
3044 }
3045
3046 gxstyl (1);
3047 gxwide (1);
3048 gas2d (pcm, pgru, 1); /* Set up scaling, draw map if apprprt */
3049
3050 gafram (pcm);
3051 idiv = 1.0; jdiv = 1.0;
3052
3053 if (flag) {
3054 gamnmx (pgrc);
3055 if (pgrc->rmin==pgrc->undef) {
3056 gaprnt (0,"Connot color vectors/barbs -- Color grid all undefined\n");
3057 flag = 0;
3058 } else gaselc(pcm,pgrc->rmin,pgrc->rmax);
3059 }
3060
3061 gamnmx (pgru);
3062 gamnmx (pgrv);
3063 if (pgru->rmin==pgru->undef) {
3064 gaprnt (0,"Cannot draw vectors/barbs -- U field all undefined\n");
3065 return;
3066 }
3067 if (pgrv->rmin==pgrv->undef) {
3068 gaprnt (0,"Cannot draw vectors/barbs -- V field all undefined\n");
3069 return;
3070 }
3071
3072 if (pcm->ccolor>=0) lcol = pcm->ccolor;
3073 else lcol=1;
3074 gxcolr (lcol);
3075 gxstyl(1);
3076 gxwide(pcm->cthick);
3077 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
3078 if (!pcm->arrflg) {
3079 umax = pgru->rmax;
3080 if (umax<fabs(pgru->rmin)) umax = fabs(pgru->rmin);
3081 vmax = pgrv->rmax;
3082 if (vmax<fabs(pgrv->rmin)) vmax = fabs(pgrv->rmin);
3083 vscal = hypot(umax,vmax);
3084 if (vscal>0.0) {
3085 x = floor(log10(vscal));
3086 y = floor(vscal/pow(10.0,x));
3087 vscal = y * pow(10.0,x);
3088 pcm->arrsiz = 0.5;
3089 } else vscal=1.0;
3090 pcm->arrmag = vscal;
3091 } else {
3092 vscal = pcm->arrmag;
3093 }
3094 pcm->arrflg = 1;
3095
3096 u = pgru->grid;
3097 v = pgrv->grid;
3098 if (flag) c = pgrc->grid;
3099
3100 /* hflg=0 idim is lat
3101 hflg=1 jdim is lat
3102 hflg=2 plot nhem
3103 hflg=3 plot shem */
3104
3105 if (pcm->hemflg==0) hflg = 2;
3106 else if (pcm->hemflg==1) hflg = 3;
3107 else {
3108 if (pgru->idim==1) hflg = 0;
3109 else if (pgru->jdim==1) hflg = 1;
3110 else {
3111 if (pcm->dmin[1] < 0.0) hflg = 3;
3112 else hflg = 2;
3113 }
3114 }
3115 for (j=1; j<=pgru->jsiz; j++) {
3116 for (i=1; i<=pgru->isiz; i++) {
3117 if (*u!=pgru->undef && *v!=pgrv->undef) {
3118 gxconv ((float)i,(float)j,&x,&y,3);
3119 gxgrmp ((float)i,(float)j,&lon,&lat);
3120 adj = gxaarw (lon, lat);
3121 hemflg = 0;
3122 if (hflg==0 && lon<0.0) hemflg = 1;
3123 if (hflg==1 && lat<0.0) hemflg = 1;
3124 if (hflg==3) hemflg = 1;
3125 if (flag) {
3126 if (*c==pgrc->undef) gxcolr(15);
3127 else {
3128 lcol = gashdc(pcm,*c);
3129 if (lcol>-1) gxcolr(lcol);
3130 }
3131 }
3132 if (lcol>-1) {
3133 if (*v==0.0 && *u==0.0) dir = 0.0;
3134 else dir = atan2(*v,*u);
3135 if (brbflg) {
3136 gabarb (x, y, pcm->digsiz*3.5, pcm->digsiz*2.0,
3137 pcm->digsiz*0.25, dir+adj, hypot(*u,*v), hemflg);
3138 } else {
3139 if (vscal>0.0) {
3140 gaarrw (x, y, dir+adj, pcm->arrsiz*hypot(*u,*v)/vscal, pcm->ahdsiz);
3141 } else {
3142 gaarrw (x, y, dir+adj, pcm->arrsiz, pcm->ahdsiz);
3143 }
3144 }
3145 }
3146 }
3147 u++; v++;
3148 if (flag) c++;
3149 }}
3150
3151 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
3152
3153 if (pcm->arlflg && vscal>0.0 && !brbflg) {
3154 gxcolr (pcm->anncol);
3155 gxwide (pcm->annthk-2);
3156 gaarrw (pcm->xsiz2-2.0,pcm->ysiz1-0.5,0.0,pcm->arrsiz, pcm->ahdsiz);
3157 sprintf (pout,"%g",vscal);
3158 len = strlen(pout);
3159 x = pcm->xsiz2 - 2.0 + (pcm->arrsiz/2.0) - 0.5*0.13*(float)len;
3160 gxchpl (pout,len,x,pcm->ysiz1-0.7,0.13,0.13,0.0);
3161 }
3162 gxwide (pcm->annthk-3);
3163 gxcolr (pcm->anncol);
3164 if (pcm->pass==0 && pcm->grdsflg)
3165 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
3166 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
3167 gaaxpl(pcm,pgru->idim,pgru->jdim);
3168 if (brbflg) gagsav (15,pcm,pgru);
3169 else gagsav (3,pcm,pgru);
3170 }
3171
3172 /* Do (for now 2-D) scatter plot */
3173
3174 int scatcol[6] = {1,3,2,4,7,8};
3175 int scattyp[6] = {1,6,4,8,7,2};
3176
gascat(struct gacmn * pcm)3177 void gascat (struct gacmn *pcm) {
3178 struct gagrid *pgr1, *pgr2;
3179 float *r1, *r2;
3180 float cmin1,cmax1,cmin2,cmax2,cint1,cint2,x,y;
3181 int siz,i,pass,im;
3182
3183 if (pcm->numgrd<2) {
3184 gaprnt (0,"Error plotting scatter field: Only one grid provided\n");
3185 return;
3186 }
3187 if (pcm->type[0]==0 || pcm->type[1]==0) {
3188 gaprnt (0,"Error plotting scatter field: stn argument(s) used\n");
3189 return;
3190 }
3191 pgr1 = pcm->result[0].pgr;
3192 pgr2 = pcm->result[1].pgr;
3193 if (pgr1->idim!=pgr2->idim || pgr1->jdim!=pgr2->jdim ||
3194 gagchk(pgr1,pgr2,pgr1->idim) ||
3195 (pgr1->jdim>-1 && gagchk(pgr1,pgr2,pgr1->jdim))) {
3196 gaprnt (0,"Error plotting scatter plot: Invalid grids\n");
3197 gaprnt (0," The two grids have difference scaling\n");
3198 return;
3199 }
3200
3201 pcm->xdim = 4;
3202 pcm->ydim = 4;
3203
3204 gamnmx (pgr1);
3205 gamnmx (pgr2);
3206 if (pgr1->rmin==pgr1->undef) {
3207 gaprnt (0,"Cannot draw scatter plot -- 1st field all undefined\n");
3208 return;
3209 }
3210 if (pgr2->rmin==pgr2->undef) {
3211 gaprnt (0,"Cannot draw scatter plot -- 2nd field all undefined\n");
3212 return;
3213 }
3214
3215 if (pcm->paflg) {
3216 pcm->xsiz1 = pcm->pxmin;
3217 pcm->xsiz2 = pcm->pxmax;
3218 pcm->ysiz1 = pcm->pymin;
3219 pcm->ysiz2 = pcm->pymax;
3220 } else {
3221 pcm->xsiz1 = 2.0;
3222 pcm->xsiz2 = pcm->xsiz - 1.5;
3223 pcm->ysiz1 = 1.00;
3224 pcm->ysiz2 = pcm->ysiz - 1.00;
3225 }
3226 gafram (pcm);
3227 gxwide (pcm->cthick);
3228 idiv = 1.0; jdiv = 1.0;
3229
3230 if (pcm->aflag != 0) {
3231 cmin1 = pcm->rmin;
3232 cmax1 = pcm->rmax;
3233 } else {
3234 cint1 = 0.0;
3235 gacsel (pgr1->rmin,pgr1->rmax,&cint1,&cmin1,&cmax1);
3236 if (cint1==0.0) {
3237 cmin1 = pgr1->rmin-5.0;
3238 cmax1 = cmin1+10.0;
3239 cint1 = 2.0;
3240 } else {
3241 cmin1 = cmin1 - 2.0*cint1;
3242 cmax1 = cmax1 + 2.0*cint1;
3243 }
3244 }
3245 if (pcm->aflag2 != 0) {
3246 cmin2 = pcm->rmin2;
3247 cmax2 = pcm->rmax2;
3248 } else {
3249 cint2 = 0.0;
3250 gacsel (pgr2->rmin,pgr2->rmax,&cint2,&cmin2,&cmax2);
3251 if (cint2==0.0) {
3252 cmin2 = pgr2->rmin-5.0;
3253 cmax2 = cmin2+10.0;
3254 cint2 = 2.0;
3255 } else {
3256 cmin2 = cmin2 - 2.0*cint2;
3257 cmax2 = cmax2 + 2.0*cint2;
3258 }
3259 }
3260
3261 printf ("%g %g %g %g \n",cmin1,cmax1,cmin2,cmax2);
3262 gxscal (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2,
3263 cmin1,cmax1,cmin2,cmax2);
3264
3265 pass = pcm->gpass[3];
3266 if (pass>5) pass=5;
3267 if (pcm->ccolor<0) gxcolr(1);
3268 else gxcolr(pcm->ccolor);
3269 gxwide (pcm->cthick);
3270 siz = pgr1->isiz*pgr1->jsiz;
3271 r1 = pgr1->grid;
3272 r2 = pgr2->grid;
3273 for (i=0; i<siz; i++) {
3274 if (*r1!=pgr1->undef && *r2!=pgr2->undef) {
3275 gxconv (*r1,*r2,&x,&y,1);
3276
3277 /*mf 950925 unable cmark for scatter plots mf*/
3278
3279 im = pcm->cmark;
3280 /* gxmark (scattyp[pass],x,y,pcm->digsiz*0.5); */
3281 gxmark (im,x,y,pcm->digsiz*0.5);
3282 }
3283 r1++; r2++;
3284 }
3285
3286 pcm->rmin = cmin1;
3287 pcm->rmax = cmax1;
3288 pcm->rint = cint1;
3289 gaaxis(1,pcm,4);
3290 pcm->rmin = cmin2;
3291 pcm->rmax = cmax2;
3292 pcm->rint = cint2;
3293 gaaxis(0,pcm,4);
3294
3295 pcm->rmin = cmin1;
3296 pcm->rmax = cmax1;
3297 pcm->aflag = 1;
3298 pcm->rmin2 = cmin2;
3299 pcm->rmax2 = cmax2;
3300 pcm->aflag2 = 1;
3301
3302 if (cmin1<0.0 && cmax1>0.0) {
3303 gxstyl(1);
3304 gxwide(3);
3305 gxconv (0.0, cmin2, &x, &y, 1);
3306 gxplot (x,y,3);
3307 gxconv (0.0, cmax2, &x, &y, 1);
3308 gxplot (x,y,2);
3309 }
3310 if (cmin2<0.0 && cmax2>0.0) {
3311 gxstyl(1);
3312 gxwide(3);
3313 gxconv (cmin1, 0.0, &x, &y, 1);
3314 gxplot (x,y,3);
3315 gxconv (cmax1, 0.0, &x, &y, 1);
3316 gxplot (x,y,2);
3317 }
3318 if (pcm->pass==0 && pcm->grdsflg)
3319 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
3320 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
3321 gagsav (3,pcm,pgr1);
3322 pcm->gpass[3]++;
3323 }
3324
3325 static float a150 = 150.0*3.1416/180.0;
3326
gaarrw(float x,float y,float ang,float siz,float asiz)3327 void gaarrw (float x, float y, float ang, float siz, float asiz) {
3328 float xx,yy;
3329
3330 if (siz<0.0001) {
3331 gxmark (2, x, y, 0.01);
3332 return;
3333 }
3334 gxplot (x,y,3);
3335 xx = x+siz*cos(ang);
3336 yy = y+siz*sin(ang);
3337 gxplot (xx,yy,2);
3338 if (asiz==0.0) return;
3339 if (asiz<0.0) asiz = -1.0*asiz*siz;
3340 gxplot (xx+asiz*cos(ang+a150),yy+asiz*sin(ang+a150),2);
3341 gxplot (xx,yy,3);
3342 gxplot (xx+asiz*cos(ang-a150),yy+asiz*sin(ang-a150),2);
3343 }
3344
3345 /* Do 2-D grid value plot */
3346
gaplvl(struct gacmn * pcm)3347 void gaplvl (struct gacmn *pcm) {
3348 struct gagrid *pgr,*pgrm;
3349 int i,j,len,lcol;
3350 float xlo,ylo,xhi,yhi,dum,*r,*m,cwid;
3351 char lab[20];
3352 int flag;
3353
3354 pgr = pcm->result[0].pgr;
3355 flag = 0;
3356 if (pcm->numgrd>1) {
3357 flag = 1;
3358 pgrm = pcm->result[1].pgr;
3359 if (pgrm->idim!=pgr->idim || pgrm->jdim!=pgr->jdim ||
3360 gagchk(pgrm,pgr,pgr->idim) || gagchk(pgrm,pgr,pgr->jdim)) {
3361 flag = 0;
3362 gaprnt (0,"Error plotting grid values: Invalid Mask grid");
3363 gaprnt (0," Mask grid ignored -- has different scaling");
3364 }
3365 }
3366
3367 if ( (pcm->rotate && (pgr->idim!=2 || pgr->jdim!=3)) ||
3368 (!pcm->rotate && pgr->idim==2 && pgr->jdim==3)) {
3369 pgr = gaflip(pgr,pcm);
3370 if (flag) pgrm = gaflip(pgrm,pcm);
3371 }
3372
3373 gxstyl (1);
3374 gxwide (1);
3375 gas2d (pcm, pgr, 1); /* Draw map and set up scaling */
3376
3377 gafram (pcm);
3378 gxwide (pcm->cthick);
3379 idiv = 1.0; jdiv = 1.0;
3380 if (pcm->ccolor>=0) lcol = pcm->ccolor;
3381 else lcol=1;
3382 gxcolr(lcol);
3383 gxstyl(1);
3384 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
3385
3386 /* Draw grid lines */
3387
3388 if (pcm->gridln != -1) {
3389 if (pcm->gridln > -1) gxcolr(pcm->gridln);
3390 for (i=1; i<=pgr->isiz; i++) {
3391 for (j=1; j<=pgr->jsiz; j++) {
3392 gxconv ((float)(i)+0.5,(float)(j)-0.5,&xlo,&ylo,3);
3393 gxconv ((float)(i)+0.5,(float)(j)+0.5,&xhi,&yhi,3);
3394 gxplot (xlo,ylo,3);
3395 gxplot (xhi,yhi,2);
3396 }
3397 }
3398
3399 for (j=1; j<=pgr->jsiz; j++) {
3400 for (i=1; i<=pgr->isiz; i++) {
3401 gxconv ((float)(i)-0.5,(float)(j)+0.5,&xlo,&ylo,3);
3402 gxconv ((float)(i)+0.5,(float)(j)+0.5,&xhi,&yhi,3);
3403 gxplot (xlo,ylo,3);
3404 gxplot (xhi,yhi,2);
3405 }
3406 }
3407 }
3408
3409 r = pgr->grid;
3410 if (flag) m = pgrm->grid;
3411 for (j=1; j<=pgr->jsiz; j++) {
3412 for (i=1; i<=pgr->isiz; i++) {
3413 if (*r!=pgr->undef) {
3414 if (flag && *m<=0.0) {
3415 gxwide (1);
3416 gxcolr (15);
3417 } else {
3418 gxwide (pcm->cthick);
3419 gxcolr (lcol);
3420 }
3421 gxconv ((float)i,(float)j,&xlo,&ylo,3);
3422 sprintf (lab,"%.*f",pcm->dignum,*r);
3423 len = strlen(lab);
3424 cwid = pcm->digsiz*(float)len;
3425 gxchln (lab,len,pcm->digsiz,&cwid);
3426 gxchpl (lab,len,xlo-cwid*0.5,ylo-pcm->digsiz*0.5,
3427 pcm->digsiz,pcm->digsiz,0.0);
3428 }
3429 r++;
3430 if (flag) m++;
3431 }}
3432
3433 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
3434 gxcolr(pcm->anncol);
3435 gxwide(pcm->annthk-3);
3436 if (pcm->pass==0 && pcm->grdsflg)
3437 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
3438 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
3439 gaaxpl(pcm,pgr->idim,pgr->jdim);
3440 gagsav (4,pcm,pgr);
3441
3442 }
3443
3444 /* Write grid to a file. */
3445
gafwrt(struct gacmn * pcm)3446 void gafwrt (struct gacmn *pcm) {
3447 struct gagrid *pgr;
3448 int size, exsz, rdw, ret, diff;
3449 int row, goff, xoff, yoff, xsiz, ysiz, incr;
3450 size_t write,written;
3451 char fmode[3];
3452
3453 pgr = pcm->result[0].pgr;
3454 size = pgr->isiz * pgr->jsiz;
3455 if (pcm->ffile == NULL ) {
3456 /* Check if ffile file exists and setting write mode */
3457 pcm->ffile = fopen(pcm->fwname,"r");
3458 if(pcm->fwappend) {
3459 strcpy(fmode,"ab");
3460 if (pcm->ffile != NULL){
3461 sprintf(pout,"Appending data to file %s.\n",pcm->fwname);
3462 gaprnt (2,pout);
3463 }
3464 } else {
3465 strcpy(fmode,"wb");
3466 if (pcm->ffile != NULL){
3467 sprintf(pout,"Replacing file %s.\n",pcm->fwname);
3468 gaprnt (2,pout);
3469 }
3470 }
3471 if (pcm->ffile != NULL) fclose(pcm->ffile);
3472 if (pcm->fwname) pcm->ffile = fopen(pcm->fwname,fmode);
3473 else pcm->ffile = fopen ("grads.fwrite",fmode);
3474 if (pcm->ffile==NULL) {
3475 gaprnt (0,"Error opening output file for fwrite\n");
3476 if (pcm->fwname) {
3477 gaprnt (0," File name is: ");
3478 gaprnt (0,pcm->fwname);
3479 gaprnt (0,"\n");
3480 } else {
3481 gaprnt (0," File name is: grads.fwrite\n");
3482 }
3483 return;
3484 }
3485 }
3486
3487 /* swap if needed. assumes 4 byte values */
3488
3489 rdw = size*4;
3490 if (BYTEORDER != pcm->fwenflg) {
3491 gabswp(pgr->grid,size);
3492 gabswp(&rdw,1);
3493 }
3494
3495 /* Handle -ex flag -- try to use exact grid coords */
3496 written = 0;
3497 exsz = size;
3498 diff = 0;
3499 /* only X or Y can vary for new fwrite code */
3500 if (pcm->fwexflg && pgr->idim<2 && pgr->jdim<2) {
3501 if (pcm->xexflg) {
3502 if (pcm->x1ex != pgr->dimmin[0]) diff=1;
3503 if (pcm->x2ex != pgr->dimmax[0]) diff=1;
3504 }
3505 if (pcm->yexflg) {
3506 if (pcm->y1ex != pgr->dimmin[1]) diff=1;
3507 if (pcm->y2ex != pgr->dimmax[1]) diff=1;
3508 }
3509 }
3510
3511 if (diff) {
3512 if (pgr->idim==0) { /* x is varying */
3513 if (pcm->xexflg) {
3514 xoff = pcm->x1ex - pgr->dimmin[0];
3515 xsiz = 1 + pcm->x2ex - pcm->x1ex;
3516 } else {
3517 xoff = 0;
3518 xsiz = pgr->isiz;
3519 }
3520 if (pgr->jdim==1 && pcm->yexflg) { /* both x and y vary */
3521 yoff = pcm->y1ex - pgr->dimmin[1];
3522 ysiz = 1 + pcm->y2ex - pcm->y1ex;
3523 } else {
3524 yoff = 0;
3525 ysiz = pgr->jsiz;
3526 }
3527 }
3528 incr = pgr->isiz;
3529 if (pgr->idim==1) { /* x is fixed; y is varying */
3530 if (pcm->yexflg) {
3531 yoff = pcm->y1ex - pgr->dimmin[1];
3532 ysiz = 1 + pcm->y2ex - pcm->y1ex;
3533 } else {
3534 yoff = 0;
3535 ysiz = pgr->isiz;
3536 }
3537 xoff = 0; xsiz = 1;
3538 incr = 1;
3539 }
3540 exsz = xsiz * ysiz;
3541 rdw = exsz*4;
3542 /* Swap the record header if necessary. fix by LIS @ NASA, 3/8/2004 ***/
3543 if (BYTEORDER != pcm->fwenflg) gabswp(&rdw,1);
3544
3545 if (pcm->fwsqflg) fwrite(&rdw,sizeof(int),1,pcm->ffile);
3546 if (pgr->idim==1) {
3547 goff = yoff;
3548 } else {
3549 goff = xoff + yoff*pgr->isiz;
3550 }
3551 for (row=0; row<ysiz; row++) {
3552 write = fwrite (pgr->grid+goff, sizeof(float), xsiz, pcm->ffile);
3553 ret = ferror(pcm->ffile);
3554 if (ret || (write != xsiz)) {
3555 sprintf(pout, "Error writing data for fwrite: %s\n", strerror(errno) );
3556 gaprnt(0, pout);
3557 }
3558 written = written + write;
3559 goff = goff + incr;
3560 }
3561 if (pcm->fwsqflg) fwrite(&rdw,sizeof(int),1,pcm->ffile);
3562 } else {
3563 if (pcm->fwsqflg) fwrite(&rdw,sizeof(int),1,pcm->ffile);
3564
3565 written = fwrite (pgr->grid, sizeof(float), size, pcm->ffile);
3566 ret = ferror(pcm->ffile);
3567 if (ret || (written != size)) {
3568 sprintf(pout, "Error writing data for fwrite: %s\n", strerror(errno) );
3569 gaprnt(0, pout);
3570 }
3571
3572 if (pcm->fwsqflg) fwrite(&rdw,sizeof(int),1,pcm->ffile);
3573 }
3574
3575 /* swap back (we don't want to change the data in case
3576 someone uses it later in the code */
3577
3578 if (BYTEORDER != pcm->fwenflg) gabswp(pgr->grid,size);
3579
3580 if (pcm->fwname) {
3581 sprintf (pout,"Wrote %i of %i elements to ", written, exsz);
3582 gaprnt (2,pout);
3583 gaprnt (2,pcm->fwname);
3584 } else {
3585 sprintf (pout,"Wrote %i of %i elements to grads.fwrite", written, exsz);
3586 gaprnt (2,pout);
3587 }
3588
3589 if (pcm->fwsqflg) gaprnt(2," as Sequential");
3590 else gaprnt(2," as Stream");
3591 if (pcm->fwenflg) gaprnt(2," Big_Endian\n");
3592 else gaprnt(2," Little_Endian\n");
3593 gagsav (20,pcm,NULL);
3594 }
3595
3596 /* Do 2-D grid fill plots */
3597
gafgrd(struct gacmn * pcm)3598 void gafgrd (struct gacmn *pcm) {
3599 struct gagrid *pgr;
3600 int i,j,k,iv,col,scol,flag,isav,ii,siz;
3601 float *xybuf,*xy,*r;
3602
3603 pgr = pcm->result[0].pgr;
3604
3605 if ( (pcm->rotate && (pgr->idim!=2 || pgr->jdim!=3)) ||
3606 (!pcm->rotate && pgr->idim==2 && pgr->jdim==3)) pgr = gaflip(pgr,pcm);
3607
3608 gas2d (pcm, pgr, 0); /* Set up scaling */
3609 idiv = 1.0; jdiv = 1.0;
3610 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
3611
3612 /* Allocate point buffer */
3613
3614 siz = (pgr->isiz+2)*4;
3615 xybuf = (float *)malloc(sizeof(float)*siz);
3616 if (xybuf==NULL) {
3617 gaprnt(0,"Memory allocation error in FGRID\n");
3618 return;
3619 }
3620 *(xybuf+siz-1) = -1;
3621
3622 /* Fill grid "boxes" */
3623
3624 r = pgr->grid;
3625 for (j=1; j<=pgr->jsiz; j++) {
3626 col = -1; scol = -1;
3627 isav = 1;
3628 for (i=1; i<=pgr->isiz; i++) {
3629 col = -1;
3630 if (*r != pgr->undef) {
3631 iv = floor(*r+0.5);
3632 for (k=0; k<pcm->fgcnt; k++) {
3633 if (iv==pcm->fgvals[k]) col = pcm->fgcols[k];
3634 }
3635 }
3636 if (col!=scol) {
3637 if (scol>-1) {
3638 xy = xybuf;
3639 for (ii=isav; ii<=i; ii++) {
3640 gxconv ((float)(ii)-0.5,(float)(j)+0.5,xy,xy+1,3);
3641 xy+=2;
3642 }
3643 for (ii=i; ii>=isav; ii--) {
3644 gxconv ((float)(ii)-0.5,(float)(j)-0.5,xy,xy+1,3);
3645 xy+=2;
3646 }
3647 *xy = *xybuf; *(xy+1) = *(xybuf+1);
3648 gxcolr(scol);
3649 gxfill (xybuf,(1+i-isav)*2+1);
3650 }
3651 isav = i;
3652 scol = col;
3653 }
3654 r++;
3655 }
3656 if (scol>-1) {
3657 xy = xybuf;
3658 for (ii=isav; ii<=pgr->isiz+1; ii++) {
3659 gxconv ((float)(ii)-0.5,(float)(j)+0.5,xy,xy+1,3);
3660 xy+=2;
3661 }
3662 for (ii=pgr->isiz+1; ii>=isav; ii--) {
3663 gxconv ((float)(ii)-0.5,(float)(j)-0.5,xy,xy+1,3);
3664 xy+=2;
3665 }
3666 *xy = *xybuf; *(xy+1) = *(xybuf+1);
3667 gxcolr(scol);
3668 gxfill (xybuf,(2+pgr->isiz-isav)*2+1);
3669 }
3670 }
3671 if (*(xybuf+siz-1) != -1) {
3672 printf ("Logic Error 16 in gafgrd. Please report error.\n");
3673 }
3674 free (xybuf);
3675 if (pgr->idim==0 && pgr->jdim==1) gawmap (pcm, 0);
3676 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
3677 gxcolr(pcm->anncol);
3678 gxwide(pcm->annthk-3);
3679 if (pcm->pass==0 && pcm->grdsflg)
3680 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
3681 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
3682 gaaxpl(pcm,pgr->idim,pgr->jdim);
3683 gafram (pcm);
3684 gagsav (5,pcm,pgr);
3685 }
3686
3687 /* Do 2-D contour plots */
3688
gacntr(struct gacmn * pcm,int filflg)3689 void gacntr (struct gacmn *pcm, int filflg) {
3690 struct gagrid *pgr;
3691 float cmin,cmax,cint,rl,rr,rrb,rmin,rmax;
3692 float *rrr,*r;
3693 int i,iexp,jexp,isz,jsz,ii,ii2,isav;
3694 char chlab[50];
3695 float pmin,pmax,ci,ccnt,cnt,tt;
3696 float slevs[100],smin;
3697 int scols[100], lcnt, smooth, irb;
3698
3699 pgr = pcm->result[0].pgr;
3700
3701 if ( (pcm->rotate && (pgr->idim!=2 || pgr->jdim!=3)) ||
3702 (!pcm->rotate && pgr->idim==2 && pgr->jdim==3)) pgr = gaflip(pgr,pcm);
3703
3704 gxstyl (1);
3705 gxwide (1);
3706 if (filflg) gas2d (pcm, pgr, 0); /* No map yet if shaded cntrs */
3707 else gas2d (pcm, pgr, 1); /* Scaling plus map */
3708
3709 /* Get contour interval */
3710
3711 gamnmx (pgr);
3712 if (pgr->rmin==pgr->undef) {
3713 gaprnt (1,"Cannot contour grid - all undefined values \n");
3714 gxcolr(1);
3715 if (pcm->dwrnflg) gxchpl ("Entire Grid Undefined",21,3.0,4.5,0.3,0.25,0.0);
3716 return;
3717 }
3718 if (!pcm->cflag) {
3719 gacsel (pgr->rmin,pgr->rmax,&(pcm->cint),&cmin,&cmax);
3720 cint = pcm->cint;
3721 if (cint==0.0) {
3722
3723 if (pcm->dwrnflg) {
3724 /*mf--- 960722
3725 new treatment of constant field -- use fgrid and red to display the grid
3726 ---mf*/
3727
3728 isav=pcm->gout2a;
3729 pcm->fgvals[0]=pgr->rmin;
3730 pcm->fgcols[0]=2;
3731 pcm->gout2a = 6;
3732 pcm->fgcnt = 1;
3733 gaplot (pcm);
3734 pcm->gout2a = isav;
3735 pcm->fgcnt = 0;
3736
3737 /*mf--- 960722 ---mf*/
3738
3739 sprintf (pout,"Constant field. Value = %g\n",pgr->rmin);
3740 i = strlen(pout);
3741 gxchpl (pout,i-1,2.0,4.5,0.27,0.23,0.0);
3742 gaprnt (1,pout);
3743 } else {
3744 sprintf (pout,"Constant field. Value = %g\n",pgr->rmin);
3745 gaprnt (1,pout);
3746 }
3747 return;
3748 }
3749 pmin = cmin;
3750 if (pmin<pcm->cmin) pmin = pcm->cmin;
3751 pmax = cmax;
3752 if (pmax>pcm->cmax) pmax = pcm->cmax;
3753 if ((pmax-pmin)/cint>100.0) {
3754 gaprnt (0,"Too many contour levels -- adjusting cint\n");
3755 while ((pmax-pmin)/cint>100.0) cint*=10.0;
3756 pcm->cint = cint;
3757 gacsel (pgr->rmin,pgr->rmax,&cint,&cmin,&cmax);
3758 }
3759 }
3760
3761 if (pcm->ccolor>=0) gxcolr(pcm->ccolor);
3762 if (pcm->ccolor<0 && pcm->rainmn==0.0 && pcm->rainmx==0.0 &&
3763 !pcm->cflag) {
3764 pcm->rainmn = cmin;
3765 pcm->rainmx = cmax;
3766 }
3767
3768 /* Expand grid if smoothing was requested */
3769
3770 idiv = 1.0; jdiv = 1.0;
3771 smooth = 0;
3772 rmin = pgr->rmin; rmax = pgr->rmax;
3773 if (pcm->csmth && (pgr->isiz<51 || pgr->jsiz<51)) {
3774 smooth = 1;
3775 iexp = 100 / pgr->isiz;
3776 jexp = 100 / pgr->jsiz;
3777 if (iexp>5) iexp = 4;
3778 if (jexp>5) jexp = 4;
3779 if (iexp<1) iexp = 1;
3780 if (jexp<1) jexp = 1;
3781 isz = ((pgr->isiz-1)*iexp) + 1;
3782 jsz = ((pgr->jsiz-1)*jexp) + 1;
3783 rrr = (float *)malloc(isz*jsz*sizeof(float));
3784 if (rrr==NULL) {
3785 gaprnt (0,"Memory Allocation Error: CSMOOTH operation\n");
3786 return;
3787 }
3788 idiv = (float)iexp;
3789 jdiv = (float)jexp;
3790 if (pcm->csmth==1) {
3791 gagexp (pgr->grid, pgr->isiz, pgr->jsiz, rrr, iexp, jexp,
3792 pgr->undef);
3793 } else {
3794 gaglin (pgr->grid, pgr->isiz, pgr->jsiz, rrr, iexp, jexp,
3795 pgr->undef);
3796 }
3797
3798 /* We may have created new contour levels. Adjust cmin and
3799 cmax appropriately */
3800
3801 if (!pcm->cflag) {
3802 rmin = 9.99E35;
3803 rmax = -9.99E35;
3804 r = rrr;
3805 cnt=0;
3806 for (i=0;i<isz*jsz;i++) {
3807 if (*r != pgr->undef) {
3808 cnt++;
3809 if (rmin>*r) rmin = *r;
3810 if (rmax<*r) rmax = *r;
3811 }
3812 r++;
3813 }
3814 if (cnt==0 || rmin==9.99e35 || rmax==-9.99e35) {
3815 rmin = pgr->undef;
3816 rmax = pgr->undef;
3817 }
3818
3819 if (rmin==pgr->undef) {
3820 gaprnt (1,"Cannot contour grid - all undefined values \n");
3821 if (pcm->dwrnflg) gxchpl ("Entire Grid Undefined",21,3.0,4.5,0.3,0.25,0.0);
3822 free(rrr);
3823 return;
3824 }
3825 while (rmin<cmin) cmin -= cint;
3826 while (rmax>cmax) cmax += cint;
3827 }
3828 }
3829 if (pcm->cflag) {
3830 gaprnt (2,"Contouring at clevs = ");
3831 for (i=0; i<pcm->cflag; i++) {
3832 sprintf (pout," %g",pcm->clevs[i]);
3833 gaprnt (2,pout);
3834 }
3835 gaprnt (2,"\n");
3836 } else {
3837 sprintf (pout,"Contouring: %g to %g interval %g \n",cmin,cmax,cint);
3838 gaprnt (2,pout);
3839 }
3840 gxclip (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2);
3841 if (pcm->rbflg) irb = pcm->rbflg - 1;
3842 else irb = 12;
3843 rrb = irb+1;
3844 if (filflg) {
3845 gaselc (pcm,pgr->rmin,pgr->rmax);
3846 if (smooth) {
3847 if (filflg==1) gxshad (rrr,isz,jsz,pcm->shdlvs,pcm->shdcls,
3848 pcm->shdcnt,pgr->undef);
3849 else gagfil (rrr,isz,jsz,pcm->shdlvs,pcm->shdcls,pcm->shdcnt,pgr->undef);
3850 } else {
3851 if (filflg==1) {
3852 gxshad (pgr->grid,pgr->isiz,pgr->jsiz,pcm->shdlvs,pcm->shdcls,
3853 pcm->shdcnt,pgr->undef);
3854 } else {
3855 gagfil (pgr->grid,pgr->isiz,pgr->jsiz,pcm->shdlvs,pcm->shdcls,
3856 pcm->shdcnt,pgr->undef);
3857 }
3858 }
3859 if (pgr->idim==0 && pgr->jdim==1) gawmap (pcm, 0);
3860 } else {
3861 gxwide (pcm->cthick);
3862 if (pcm->cflag) {
3863 for (i=0; i<pcm->cflag; i++) {
3864 rr = pcm->clevs[i];
3865 if (rr<0.0&&pcm->cstyle==-9) gxstyl (3);
3866 else gxstyl(pcm->cstyle);
3867 if (pcm->ccolor < 0 && pcm->ccflg == 0) {
3868 if (pcm->cflag==1) ii=irb/2;
3869 else ii = (int)((float)(i*irb)/((float)(pcm->cflag-1)));
3870 if (ii>irb) ii=irb;
3871 if (ii<0) ii=0;
3872 if (pcm->rbflg>0) {
3873 if (pcm->ccolor==-1) gxcolr(pcm->rbcols[ii]);
3874 else gxcolr(pcm->rbcols[irb-ii]);
3875 } else {
3876 if (pcm->ccolor==-1) gxcolr(rcols[ii]);
3877 else gxcolr(rcols[12-ii]);
3878 }
3879 }
3880 if (pcm->ccflg) {
3881 ii = i;
3882 if (ii>=pcm->ccflg) ii = pcm->ccflg-1;
3883 gxcolr (pcm->ccols[ii]);
3884 }
3885 if (pcm->clstr) sprintf(chlab,pcm->clstr,rr);
3886 else sprintf (chlab,"%g",rr);
3887 chlab[15] = '\0'; /* Limit length. this is bad */
3888 if (smooth) {
3889 gxclev (rrr,isz,jsz,1,isz,1,jsz,rr,pgr->undef,
3890 chlab,pcm->cterp,pcm->clab);
3891 } else {
3892 gxclev (pgr->grid,pgr->isiz,pgr->jsiz,1,pgr->isiz,1,
3893 pgr->jsiz,rr,pgr->undef,chlab,pcm->cterp,pcm->clab);
3894 }
3895 }
3896 } else {
3897 for (rl=cmin;rl<=cmax+(cint/2.0);rl+=cint) {
3898 if (rl<pcm->cmin || rl>pcm->cmax) continue;
3899 if (pcm->blkflg&&rl>=pcm->blkmin&&rl<=pcm->blkmax) continue;
3900 if (rl<0.0) ii = (int)((rl/cint)-0.1);
3901 else ii = (int)((rl/cint)+0.1);
3902 rr = (float)ii * cint;
3903 if (rr<0.0&&pcm->cstyle==-9) gxstyl (3);
3904 else gxstyl(pcm->cstyle);
3905 if (pcm->ccolor < 0) {
3906 ii = (int)(rrb*(rr-pcm->rainmn)/(pcm->rainmx-pcm->rainmn));
3907 if (ii>irb) ii=irb;
3908 if (ii<0) ii=0;
3909 if (pcm->rbflg>0) {
3910 if (pcm->ccolor==-1) gxcolr(pcm->rbcols[ii]);
3911 else gxcolr(pcm->rbcols[irb-ii]);
3912 } else {
3913 if (pcm->ccolor==-1) gxcolr(rcols[ii]);
3914 else gxcolr(rcols[12-ii]);
3915 }
3916 }
3917 tt = rl/(cint*(float)pcm->clskip);
3918
3919 /*mf add fuzz for cray mf*/
3920
3921 ii = (int)(tt+0.009);
3922 ii2 = (int)(tt-0.009);
3923 if (fabs(tt-(float)ii)<0.01 || fabs(tt-(float)ii2)<0.01 ) {
3924 if (pcm->clstr) {
3925 sprintf(chlab,pcm->clstr,rr);
3926 } else {
3927 sprintf (chlab,"%g",rr);
3928 }
3929 chlab[15] = '\0'; /* Limit length. this is bad */
3930 } else {
3931 chlab[0] = '\0';
3932 }
3933 if (smooth) {
3934 gxclev (rrr,isz,jsz,1,isz,1,jsz,rr,pgr->undef,
3935 chlab,pcm->cterp,pcm->clab);
3936 } else {
3937 gxclev (pgr->grid,pgr->isiz,pgr->jsiz,1,pgr->isiz,1,
3938 pgr->jsiz,rr,pgr->undef,chlab,pcm->cterp,pcm->clab);
3939 }
3940 }
3941 }
3942 if (pcm->clcol>-1) gxcolr (pcm->clcol);
3943 if (pcm->clthck>-1) gxwide (pcm->clthck);
3944 gxclab (pcm->clsiz,pcm->clab,pcm->clcol);
3945 gxcrel ();
3946 }
3947 if (smooth) free(rrr);
3948 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
3949 gxcolr(pcm->anncol);
3950 gxwide (pcm->annthk-3);
3951 if (pcm->pass==0 && pcm->grdsflg)
3952 gxchpl("GrADS: COLA/IGES",16,0.05,0.05,0.1,0.08,0.0);
3953 if (pcm->pass==0 && pcm->timelabflg) gatmlb(pcm);
3954 gaaxpl(pcm,pgr->idim,pgr->jdim);
3955 gafram (pcm);
3956 if (filflg==1) gagsav (2,pcm,pgr);
3957 else if (filflg==2) gagsav (16,pcm,pgr);
3958 else gagsav (1,pcm,pgr);
3959 }
3960
3961 /* Routine to perform grid level scaling. The address of this
3962 routine is passed to gxgrid. */
3963
gaconv(float s,float t,float * x,float * y)3964 void gaconv (float s, float t, float *x, float *y) {
3965 s = ((s-1.0)/idiv)+1.0;
3966 t = ((t-1.0)/jdiv)+1.0;
3967 if (iconv==NULL) *x = s+ioffset;
3968 else *x = iconv(ivars, (s+ioffset));
3969 if (jconv==NULL) *y = t+joffset;
3970 else *y = jconv(jvars, (t+joffset));
3971 }
3972
3973 /* Draw axis labels. axis = 1 is x axis, axis = 0 is y axis */
3974
gaaxis(int axis,struct gacmn * pcm,int dim)3975 void gaaxis (int axis, struct gacmn *pcm, int dim) {
3976 struct gafile *pfi;
3977 float m,b,x,y,cs,xx,cwid;
3978 float v,vmin,vmax,vincr,vstrt,vend,tt,pos;
3979 float *tvals;
3980 int ii,len,tinc,is,ic,colr,thck,flg,i,cnt;
3981 struct dt tstrt,tincr,addmo,twrk;
3982 char lab[30],olab[30],*chlb;
3983
3984 if (axis && pcm->xlab==0) return;
3985 if (!axis && pcm->ylab==0) return;
3986 addmo.yr = 0L;
3987 addmo.mo = 1L;
3988 addmo.dy = 0L;
3989 addmo.hr = 0L;
3990 addmo.mn = 0L;
3991 if (axis==1) {
3992 cs = pcm->xlsiz;
3993 colr = pcm->xlcol;
3994 thck = pcm->xlthck;
3995 if (pcm->xlside) pos = pcm->ysiz2 + pcm->xlpos;
3996 else pos = pcm->ysiz1 + pcm->xlpos;
3997 if (pcm->xlpos!=0.0) {
3998 gxcolr (pcm->anncol);
3999 gxwide (pcm->annthk);
4000 gxstyl (1);
4001 gxplot (pcm->xsiz1,pos,3);
4002 gxplot (pcm->xsiz2,pos,2);
4003 }
4004 } else {
4005 cs = pcm->ylsiz;
4006 colr = pcm->ylcol;
4007 thck = pcm->ylthck;
4008 if (pcm->ylside) pos = pcm->xsiz2 + pcm->ylpos;
4009 else pos = pcm->xsiz1 + pcm->ylpos;
4010 if (pcm->ylpos!=0.0) {
4011 gxcolr (pcm->anncol);
4012 gxwide (pcm->annthk);
4013 gxstyl (1);
4014 gxplot (pos,pcm->ysiz1,3);
4015 gxplot (pos,pcm->ysiz2,2);
4016 }
4017 }
4018 /* Select axis min and max */
4019
4020 vincr=0.0;
4021 if (dim==4) {
4022 vmin = pcm->rmin; vmax = pcm->rmax;
4023 }else if (dim==3) {
4024 pfi = pcm->pfid;
4025 tvals = pfi->abvals[3];
4026 vmin = t2gr(tvals,&(pcm->tmin));
4027 vmax = t2gr(tvals,&(pcm->tmax));
4028 } else {
4029 vmin = pcm->dmin[dim];
4030 vmax = pcm->dmax[dim];
4031 }
4032 if (axis && pcm->xlabs) {
4033 vmin = 1.0;
4034 vmax = (float)pcm->ixlabs;
4035 vincr = 1.0;
4036 dim = 5;
4037 }
4038 if (!axis && pcm->ylabs) {
4039 vmin = 1.0;
4040 vmax = (float)pcm->iylabs;
4041 vincr = 1.0;
4042 dim = 5;
4043 }
4044 if (axis && pcm->axflg && (dim!=2 || pcm->zlog==0)) {
4045 vmin = pcm->axmin;
4046 vmax = pcm->axmax;
4047 vincr = pcm->axint;
4048 dim = 5;
4049 gaprnt (1,"Warning: X axis labels overridden by SET XAXIS.\n");
4050 gaprnt (1," Labels may not reflect correct scaling for dimensions or data.\n");
4051 }
4052 if (!axis && pcm->ayflg && (dim!=2 || pcm->zlog==0)) {
4053 vmin = pcm->aymin;
4054 vmax = pcm->aymax;
4055 vincr = pcm->ayint;
4056 dim = 5;
4057 gaprnt (1,"Warning: Y axis labels overridden by SET YAXIS.\n");
4058 gaprnt (1," Labels may not reflect correct scaling for dimensions or data.\n");
4059 }
4060 if (vmin==vmax) {
4061 gaprnt (0,"gaaxis internal logic check 24 \n");
4062 return;
4063 }
4064
4065 /* Handle axis flipping */
4066
4067 if (axis) {
4068 if (pcm->xflip) {
4069 m=(pcm->xsiz2-pcm->xsiz1)/(vmin-vmax);
4070 b=pcm->xsiz1-(m*vmax);
4071 } else {
4072 m=(pcm->xsiz2-pcm->xsiz1)/(vmax-vmin);
4073 b=pcm->xsiz1-(m*vmin);
4074 }
4075 } else {
4076 if (pcm->yflip) {
4077 m=(pcm->ysiz2-pcm->ysiz1)/(vmin-vmax);
4078 b=pcm->ysiz1-(m*vmax);
4079 } else {
4080 m=(pcm->ysiz2-pcm->ysiz1)/(vmax-vmin);
4081 b=pcm->ysiz1-(m*vmin);
4082 }
4083 }
4084
4085 /* Select label interval */
4086
4087 if (vmin>vmax) {
4088 v = 1.0*vmax; /* Avoid optimization */
4089 vmax = vmin;
4090 vmin = v;
4091 }
4092 if (dim==3) {
4093 tinc = gatinc (pcm, &tstrt, &tincr);
4094 } else {
4095 flg = 1;
4096 if (axis==1 && pcm->xlint!=0.0) {vincr=pcm->xlint; flg=0;}
4097 if (axis==0 && pcm->ylint!=0.0) {vincr=pcm->ylint; flg=0;}
4098 if (vincr<0.0) {
4099 vincr = -1.0 * vincr;
4100 vstrt = vmin;
4101 vend = vmax;
4102 vend = vend+(vincr*0.5);
4103 } else {
4104 gacsel (vmin,vmax,&vincr,&vstrt,&vend);
4105 if (vincr==0.0) {
4106 gaprnt (0,"gaaxis internal logic check 25\n");
4107 return;
4108 }
4109 if (dim==1 && flg) {
4110 if (vincr>19.9) vincr=30.0;
4111 else if (vincr>10.0) vincr=10.0;
4112 }
4113 if (dim==0 && flg) {
4114 if (vincr>74.5) vincr=90.0;
4115 else if (vincr>44.5) vincr=60.0;
4116 else if (vincr>24.9) vincr=30.0;
4117 else if (vincr>14.5) vincr=20.0;
4118 else if (vincr>10.0) vincr=10.0;
4119 }
4120 gacsel (vmin,vmax,&vincr,&vstrt,&vend);
4121 vend = vend+(vincr*0.5);
4122 }
4123 }
4124
4125 gxcolr(colr);
4126 gxwide(thck);
4127 gxstyl(1);
4128 if (dim!=3) {
4129 if (axis==1 && pcm->xlflg>0) cnt = pcm->xlflg;
4130 else if (axis==0 && pcm->ylflg>0) cnt = pcm->ylflg;
4131 else {
4132 cnt = 1.0 + (vend-vstrt)/vincr;
4133 if (cnt>50) cnt=50;
4134 for (i=0; i<cnt; i++) {
4135 v = vstrt+vincr*(float)i;
4136 if (fabs(v/vincr)<1e-5) v=0.0;
4137 if (axis) *(pcm->xlevs+i) = v;
4138 else *(pcm->ylevs+i) = v;
4139 }
4140 }
4141
4142 i = 0;
4143 if (axis==1 && pcm->xlabs) chlb = pcm->xlabs;
4144 if (axis==0 && pcm->ylabs) chlb = pcm->ylabs;
4145 while (i<cnt) {
4146 if (axis) v = *(pcm->xlevs+i);
4147 else v = *(pcm->ylevs+i);
4148 if (axis==1 && pcm->xlstr) sprintf(lab,pcm->xlstr,v);
4149 else if (axis==0 && pcm->ylstr) sprintf(lab,pcm->ylstr,v);
4150 else if ( ( axis==1 && pcm->xlabs ) || ( axis==0 && pcm->ylabs) ) {
4151 sprintf(lab,chlb,v);
4152 while (*chlb) chlb++;
4153 chlb++;
4154 }
4155 else {
4156 if (dim==0 && pcm->mproj>0) len = galnch(v,lab);
4157 else if (dim==1 && pcm->mproj>0) len = galtch(v,lab);
4158 else sprintf(lab,"%g",v);
4159 }
4160 len=0;
4161 while (lab[len]) len++;
4162 cwid = (float)len*cs;
4163 gxchln (lab,len,cs,&cwid);
4164 if (axis) {
4165 x = (v*m)+b;
4166 if (dim==2 && pcm->zlog) gxconv(v,pos,&x,&tt,2);
4167 else if (dim==1 && pcm->coslat) gxconv(v,pos,&x,&tt,2);
4168 if (x<pcm->xsiz1-0.05 || x>pcm->xsiz2+0.05) {
4169 i++;
4170 continue;
4171 }
4172 gxplot (x,pos,3);
4173 if (pcm->xlside) gxplot (x,pos+(cs*0.4),2);
4174 else gxplot (x,pos-(cs*0.4),2);
4175 if (pcm->grflag==1 || pcm->grflag==3) {
4176 gxwide (1);
4177 gxstyl (pcm->grstyl);
4178 gxcolr (pcm->grcolr);
4179 gxplot (x,pcm->ysiz1,3);
4180 gxplot (x,pcm->ysiz2,2);
4181 gxcolr(colr);
4182 gxwide(thck);
4183 gxstyl(1);
4184 }
4185 x = x - cwid*0.4;
4186 if (pcm->xlside) y = pos + (cs*0.7);
4187 else y = pos - (cs*1.7);
4188 } else {
4189 y = (v*m)+b;
4190 if (dim==2 && pcm->zlog) gxconv(pcm->xsiz1,v,&tt,&y,2);
4191 else if (dim==1 && pcm->coslat) gxconv(pcm->xsiz1,v,&tt,&y,2);
4192 if (y<pcm->ysiz1-0.05 || y>pcm->ysiz2+0.05) {
4193 i++;
4194 continue;
4195 }
4196 gxplot (pos,y,3);
4197 if (pcm->ylside) {
4198 gxplot (pos+(cs*0.4),y,2);
4199 x = pos + cs*0.8;
4200 } else {
4201 gxplot (pos-(cs*0.4),y,2);
4202 x = pos - (cwid+cs)*0.8;
4203 if (pcm->yllow<(cwid+cs)*0.8) pcm->yllow = (cwid+cs)*0.8;
4204 }
4205 if (pcm->grflag==1 || pcm->grflag==2) {
4206 gxwide (1);
4207 gxstyl (pcm->grstyl);
4208 gxcolr (pcm->grcolr);
4209 gxplot (pcm->xsiz1,y,3);
4210 gxplot (pcm->xsiz2,y,2);
4211 gxwide(thck);
4212 gxcolr(colr);
4213 gxstyl(1);
4214 }
4215 y = y - (cs*0.5);
4216 }
4217 gxchpl(lab,len,x,y,cs,cs*0.8,0.0);
4218 lab[9] = '\0';
4219 i++;
4220 }
4221 } else {
4222
4223 /* Do Date/Time labeling */
4224
4225 strcpy (olab,"mmmmmmmmmmmmmmmm");
4226 while (timdif(&tstrt,&(pcm->tmax))>-1L) {
4227 len = gat2ch(&tstrt,tinc,lab);
4228 v = t2gr(tvals,&tstrt);
4229 if (axis) {
4230 x = (v*m)+b;
4231 gxplot (x,pos,3);
4232 if (pcm->xlside) gxplot (x,pos+(cs*0.4),2);
4233 else gxplot (x,pos-(cs*0.4),2);
4234
4235 /*mf 971001 this is a bug if (pcm->grflag) { mf */
4236
4237 if (pcm->grflag==1 || pcm->grflag==3) {
4238 gxwide (1);
4239 /*
4240 if (pcm->grstyl==5) gxwide(3);
4241 */
4242 gxstyl (pcm->grstyl);
4243 gxcolr (pcm->grcolr);
4244 gxplot (x,pcm->ysiz1,3);
4245 gxplot (x,pcm->ysiz2,2);
4246 gxwide(thck);
4247 gxcolr(colr);
4248 gxstyl(1);
4249 }
4250 if (pcm->xlside) y = pos + (cs*0.7);
4251 else y = pos - (cs*1.7);
4252 ii = 0;
4253 if (tinc>3) {
4254 if (tinc==5) len = 6;
4255 if (tinc==4) len = 3;
4256 if (cmpch(&(lab[ii]),&(olab[ii]),len)) {
4257 cwid = len*cs;
4258 gxchln(&lab[ii],len,cs,&cwid);
4259 xx = x - cwid*0.4;
4260 gxchpl(&lab[ii],len,xx,y,cs,cs*0.8,0.0);
4261 }
4262 if (pcm->xlside) y = y + cs*1.4;
4263 else y = y - cs*1.4;
4264 ii = len;
4265 }
4266 if (tinc>1 && pcm->tlsupp<2) {
4267 len = 5;
4268 if (tinc==2) len=3;
4269 if (cmpch(&(lab[ii]),&(olab[ii]),len)) {
4270 if (lab[ii]=='0') {ii++; len--;}
4271 cwid = len*cs;
4272 gxchln(&lab[ii],len,cs,&cwid);
4273 xx = x - cwid*0.4;
4274 gxchpl(&lab[ii],len,xx,y,cs,cs*0.8,0.0);
4275 }
4276 if (pcm->xlside) y = y + cs*1.4;
4277 else y = y - cs*1.4;
4278 ii += len;
4279 }
4280 len = 4;
4281 if (tinc!=2 || tstrt.yr!=9999) {
4282 if (pcm->tlsupp==0) {
4283 if (cmpch(&(lab[ii]),&(olab[ii]),len)) {
4284 cwid = len*cs;
4285 gxchln(&lab[ii],len,cs,&cwid);
4286 xx = x - cwid*0.4;
4287 gxchpl(&lab[ii],len,xx,y,cs,cs*0.8,0.0);
4288 }
4289 }
4290 }
4291 strcpy (olab,lab);
4292 } else {
4293 y = (v*m)+b;
4294 gxplot (pos,y,3);
4295 if (pcm->ylside) gxplot (pos+(cs*0.4),y,2);
4296 else gxplot (pos-(cs*0.4),y,2);
4297 if (pcm->grflag==1 || pcm->grflag==2) {
4298 gxwide (1);
4299 /*
4300 if (pcm->grstyl==5) gxwide(3);
4301 */
4302 gxstyl (pcm->grstyl);
4303 gxcolr (pcm->grcolr);
4304 gxplot (pcm->xsiz1,y,3);
4305 gxplot (pcm->xsiz2,y,2);
4306 gxwide(thck);
4307 gxcolr(colr);
4308 gxstyl(1);
4309 }
4310 ii = 0;
4311 if (pcm->tlsupp==1) len = len - 4;
4312 if (pcm->tlsupp==2) len = len - 9;
4313 if (len > 1) {
4314 if (tinc==3&&lab[0]=='0') {ii=1; len--;}
4315 y = y - (cs*0.5);
4316 cwid = len*cs;
4317 gxchln(&lab[ii],len,cs,&cwid);
4318 if (pcm->ylside) {
4319 x = pos + cs*0.8;
4320 } else {
4321 x = pos - (cwid+cs)*0.8;
4322 if (pcm->yllow<(cwid+cs)*0.8) pcm->yllow = (cwid+cs)*0.8;
4323 }
4324 gxchpl(&lab[ii],len,x,y,cs,cs*0.8,0.0);
4325 }
4326 }
4327
4328 /* Get next date/time. */
4329
4330 twrk = tincr;
4331 timadd (&tstrt, &twrk);
4332 tstrt = twrk;
4333 if (tincr.dy>1L&&(tstrt.dy==31L||(tstrt.dy==29L&&tstrt.dy==2))) {
4334 tstrt.dy = 1L;
4335 twrk = addmo;
4336 timadd (&tstrt, &twrk);
4337 tstrt = twrk;
4338 }
4339 if (tincr.dy>3&&tstrt.dy==3&&(tstrt.dy==2||tstrt.dy==3)) tstrt.dy = 1;
4340 }
4341 }
4342 }
4343
4344 /* Set up map projection scaling. */
4345
gamscl(struct gacmn * pcm)4346 void gamscl (struct gacmn *pcm) {
4347 float aspect,ltdif,lndif,aspect2;
4348 float xdif,ydif,xlo,xhi,ylo,yhi;
4349 int rc, flag, proj;
4350
4351 if (pcm->paflg) {
4352 mpj.xmn = pcm->pxmin;
4353 mpj.xmx = pcm->pxmax;
4354 mpj.ymn = pcm->pymin;
4355 mpj.ymx = pcm->pymax;
4356 } else {
4357 mpj.xmn = 0.5;
4358 mpj.xmx = pcm->xsiz - 0.5;
4359 mpj.ymn = 0.75;
4360 mpj.ymx = pcm->ysiz - 0.75;
4361 }
4362 if (pcm->mpflg==4 && pcm->mproj>2) {
4363 mpj.lnmn = pcm->mpvals[0]; mpj.lnmx = pcm->mpvals[1];
4364 mpj.ltmn = pcm->mpvals[2]; mpj.ltmx = pcm->mpvals[3];
4365 } else {
4366 mpj.lnmn = pcm->dmin[0]; mpj.lnmx = pcm->dmax[0];
4367 mpj.ltmn = pcm->dmin[1]; mpj.ltmx = pcm->dmax[1];
4368 }
4369 flag = 0;
4370 if (pcm->mproj==3) {
4371 rc = gxnste (&mpj);
4372 if (rc) {
4373 gaprnt (0,"Map Projection Error: Invalid coords for NPS\n");
4374 gaprnt (0," Will use linear lat-lon projection\n");
4375 flag = 1;
4376 }
4377 } else if (pcm->mproj==4) {
4378 rc = gxsste (&mpj);
4379 if (rc) {
4380 gaprnt (0,"Map Projection Error: Invalid coords for SPS\n");
4381 gaprnt (0," Will use linear lat-lon projection\n");
4382 flag = 1;
4383 }
4384 } else if (pcm->mproj==5) {
4385 rc = gxrobi (&mpj);
4386 if (rc) {
4387 gaprnt (0,"Map Projection Error: Invalid coords for Robinson\n");
4388 gaprnt (0," Will use linear lat-lon projection\n");
4389 flag = 1;
4390 }
4391
4392 /*------- DKRZ appends Projections
4393 10.08.95 Karin Meier (karin.meier@dkrz.de) -------*/
4394
4395 } else if (pcm->mproj==6) {
4396 rc = gxmoll (&mpj);
4397 if (rc) {
4398 gaprnt (0,"Map Projection Error: Invalid coords for Mollweide\n");
4399 gaprnt (0," Will use linear lat-lon projection\n");
4400 flag = 1;
4401 }
4402 } else if (pcm->mproj==7) {
4403 rc = gxortg (&mpj);
4404 if (rc) {
4405 gaprnt (0,"Map Projection Error: Invalid coords for Orthographic Projection\n");
4406 gaprnt (0," Will use linear lat-lon projection\n");
4407 flag = 1;
4408 }
4409 } else if (pcm->mproj==13) {
4410 rc = gxlamc (&mpj);
4411 if (rc) {
4412 gaprnt (0,"Map Projection Error: Invalid coords for Lambert conformal Projection\n");
4413 gaprnt (0," Will use linear lat-lon projection\n");
4414 flag = 1;
4415 }
4416 }
4417
4418
4419
4420 if (pcm->mproj==2 || flag) rc = gxltln (&mpj);
4421 else if (pcm->mproj<2) rc = gxscld (&mpj, pcm->xflip, pcm->yflip);
4422 if (rc) {
4423 gaprnt (0,"Map Projection Error: Internal Logic check\n");
4424 return;
4425 };
4426
4427 pcm->xsiz1 = mpj.axmn;
4428 pcm->xsiz2 = mpj.axmx;
4429 pcm->ysiz1 = mpj.aymn;
4430 pcm->ysiz2 = mpj.aymx;
4431 }
4432
4433 /* Plot map with proper data set, line style, etc. */
4434 /* If iflg true, check pass number */
4435
gawmap(struct gacmn * pcm,int iflg)4436 void gawmap (struct gacmn *pcm, int iflg) {
4437 struct mapopt mopt;
4438 int i;
4439
4440 if (pcm->mproj==0 || (pcm->pass > 0 && iflg)) return;
4441 if (pcm->mpdraw==0) return;
4442 gxclip (pcm->xsiz1, pcm->xsiz2, pcm->ysiz1, pcm->ysiz2);
4443 if (pcm->mapcol<0) {
4444 if (iflg) {
4445 mopt.dcol = pcm->mcolor;
4446 mopt.dstl = 1;
4447 mopt.dthk = 1;
4448 } else {
4449 mopt.dcol = 1;
4450 mopt.dstl = 1;
4451 mopt.dthk = 4;
4452 }
4453 } else {
4454 mopt.dcol = pcm->mapcol;
4455 mopt.dstl = pcm->mapstl;
4456 mopt.dthk = pcm->mapthk;
4457 }
4458 mopt.lnmin = pcm->dmin[0]; mopt.lnmax = pcm->dmax[0];
4459 mopt.ltmin = pcm->dmin[1]; mopt.ltmax = pcm->dmax[1];
4460 mopt.mcol = pcm->mpcols;
4461 mopt.mstl = pcm->mpstls;
4462 mopt.mthk = pcm->mpthks;
4463 i = 0;
4464 while (pcm->mpdset[i]) {
4465 mopt.mpdset = pcm->mpdset[i];
4466 gxdmap (&(mopt));
4467 i++;
4468 }
4469 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
4470 }
4471
4472 /* Select a contour interval. */
4473
gacsel(float rmin,float rmax,float * cint,float * cmin,float * cmax)4474 void gacsel (float rmin, float rmax, float *cint, float *cmin,
4475 float *cmax) {
4476 float rdif,norml,w1,w2,t1,t2;
4477 int power;
4478
4479 rdif = rmax-rmin;
4480 if (rdif==0.0) {
4481 *cint=0.0;
4482 *cmin=0.0;
4483 *cmax=0.0;
4484 return;
4485 }
4486
4487 if (*cint==0.0) {
4488 rdif = rdif/10.0; /* appx. 10 contour intervals */
4489 #ifndef __sgi
4490 w2 = (float)floor(log10((double)rdif));
4491 w1 = (float)pow((double)10.0,(double)w2);
4492 #else
4493 w2 = floorf(log10f(rdif));
4494 w1 = powf(10.0,w2);
4495 #endif
4496 norml = rdif/w1; /* normalized contour interval */
4497
4498 /*mf -------- mf*/
4499 if(norml < 1.0 ) norml=1.0 ; /* a 1.0=0.9999999999 problem on the cray */
4500 /*mf -------- mf*/
4501
4502 if (norml>=1.0&&norml<=1.5) *cint=1.0;
4503 else if (norml>1.5&&norml<=2.5) *cint=2.0;
4504 else if (norml>2.5&&norml<=3.5) *cint=3.0;
4505 else if (norml>3.5&&norml<=7.5) *cint=5.0;
4506 else *cint=10.0;
4507
4508 *cint = *cint*w1;
4509 }
4510
4511 *cmin = *cint * ceil(rmin/(*cint));
4512 *cmax = *cint * floor(rmax/(*cint));
4513
4514 /* Check for interval being below machine epsilon for these
4515 values */
4516
4517 t1 = rmin + *cint;
4518 t2 = rmax - *cint;
4519 if (t1 == rmin || t2 == rmax) {
4520 *cint=0.0;
4521 *cmin=0.0;
4522 *cmax=0.0;
4523 }
4524 }
4525
4526 /* Convert longitude to character */
4527
galnch(float lon,char * ch)4528 int galnch (float lon, char *ch) {
4529 int len;
4530
4531 while (lon<=-180.0) lon+=360.0;
4532 while (lon>180.0) lon-=360.0;
4533
4534 sprintf(ch,"%g",fabs(lon));
4535 len=0;
4536 while (ch[len]) len++;
4537 if (lon<0.0) {
4538 ch+=len;
4539 *ch='W';
4540 *(ch+1)='\0';
4541 len++;
4542 }
4543 if (lon>0.0&&lon<180.0) {
4544 ch+=len;
4545 *ch='E';
4546 *(ch+1)='\0';
4547 len++;
4548 }
4549 return (len);
4550 }
4551
4552 /* Convert latitude to character */
4553
galtch(float lat,char * ch)4554 int galtch (float lat, char *ch) {
4555
4556 int len;
4557
4558 sprintf(ch,"%g",fabs(lat));
4559 len=0;
4560 while (ch[len]) len++;
4561 if (lat<0.0) {
4562 ch+=len;
4563 *ch='S';
4564 *(ch+1)='\0';
4565 len++;
4566 } else if (lat>0.0) {
4567 ch+=len;
4568 *ch='N';
4569 *(ch+1)='\0';
4570 len++;
4571 } else {
4572 *ch='E';
4573 *(ch+1)='Q';
4574 *(ch+2)='\0';
4575 len=2;
4576 }
4577 return (len);
4578 }
4579
4580 /* Expand a grid using bi-linear interpolation. Same args as
4581 gagexp. */
4582
gaglin(float * g1,int cols,int rows,float * g2,int exp1,int exp2,float mval)4583 void gaglin (float *g1, int cols, int rows, float *g2,
4584 int exp1, int exp2, float mval ) {
4585 float v1,v2,xoff,yoff,x,y,xscl,yscl;
4586 int i,j,p1,p2,p3,p4,isiz,jsiz;
4587
4588 isiz = (cols-1)*exp1+1;
4589 jsiz = (rows-1)*exp2+1;
4590 printf ("isiz, jsiz = %i %i \n",isiz,jsiz);
4591 xscl = (float)(cols-1)/(float)(isiz-1);
4592 yscl = (float)(rows-1)/(float)(jsiz-1);
4593 for (j=0; j<jsiz; j++) {
4594 for (i=0; i<isiz; i++) {
4595 x = (float)i*xscl;
4596 y = (float)j*yscl;
4597 if (x<0.0) x=0.0;
4598 if (y<0.0) y=0.0;
4599 if (x>=(float)cols) x=(float)cols-0.001; /* fudge the edges */
4600 if (y>=(float)rows) y=(float)rows-0.001;
4601 p1 = (int)x + cols*(int)y;
4602 p2 = p1+1;
4603 p4 = p1+cols;
4604 p3 = p4+1;
4605 if (g1[p1]==mval || g1[p2]==mval || g1[p3]==mval || g1[p4]==mval) {
4606 *g2 = mval;
4607 } else {
4608 xoff = x-floor(x);
4609 yoff = y-floor(y);
4610 v1 = g1[p1] + (g1[p2]-g1[p1])*xoff;
4611 v2 = g1[p4] + (g1[p3]-g1[p4])*xoff;
4612 *g2 = v1 + (v2-v1)*yoff;
4613 }
4614 g2++;
4615 } }
4616 }
4617
4618 /* void gagexp (float * g1, int cols, int rows, float *g2,
4619 int exp1, int exp2, float mval);
4620
4621 Routine to expand a grid (Grid EXPand) into a larger grid
4622 using bi-directional cubic interpolation. It is expected that this
4623 routine will be used to make a finer grid for a more pleasing
4624 contouring result.
4625
4626 g1 - addr of smaller grid
4627 cols - number of columns in smaller grid
4628 rows - number of rows in smaller grid
4629 g2 - addr of where larger grid is to go
4630 exp1 - expansion factor for rows (1=no expansion, 20 max);
4631 exp2 - expansion factor for columns (1=no expansion, 20 max);
4632 mval - missing data value
4633
4634 number of output columns = (cols-1)*exp1+1
4635 number of output rows = (rows-1)*exp2+1 */
4636
4637
gagexp(float * g1,int cols,int rows,float * g2,int exp1,int exp2,float mval)4638 void gagexp (float * g1, int cols, int rows, float *g2,
4639 int exp1, int exp2, float mval ) {
4640
4641 float *p1, *p2;
4642 float s1,s2,a,b;
4643 float pows1[20],pows2[20],pows3[20];
4644 int ii,jj,k,c,d,e;
4645 float t;
4646 float kurv;
4647
4648 kurv=1.0; /* Curviness factor for spline */
4649
4650 /* First go through each row and interpolate the missing columns.
4651 Edge conditions (sides and missing value boundries) are handled
4652 by assuming a linear slope at the edge. */
4653
4654 for (k=0; k<exp1; k++) {
4655 t=( (float)k ) / ( (float)exp1 ) ;
4656 pows1[k]=t;
4657 pows2[k]=t*t;
4658 pows3[k]=t*pows2[k];
4659 }
4660
4661 p1=g1; p2=g2;
4662 c=((cols-1)*exp1)+1;
4663 d=c*exp2;
4664 e=c*(exp2-1);
4665
4666 for (jj=0; jj<rows; jj++) {
4667 for (ii=0; ii<cols-1; ii++,p1++) {
4668 if (*p1==mval || *(p1+1)==mval) {
4669 *p2 = *p1;
4670 p2++;
4671 for (k=1; k<exp1; k++,p2++) *p2=mval;
4672 } else {
4673 if (ii==0 || *(p1-1)==mval) s1 = *(p1+1) - *p1;
4674 else s1 = ( *(p1+1) - *(p1-1) )*0.5;
4675 if (ii==(cols-2) || *(p1+2)==mval) s2 = *(p1+1) - *p1;
4676 else s2 = ( *(p1+2) - *p1 )*0.5;
4677 s1*=kurv; s2*=kurv;
4678 a = s1 + 2.0*(*p1) - 2.0*(*(p1+1)) + s2;
4679 b = 3.0*(*(p1+1)) - s2 - 2.0*s1 - 3.0*(*p1);
4680 *p2 = *p1;
4681 p2++;
4682 for (k=1; k<exp1; k++,p2++) {
4683 *p2 = a*pows3[k] + b*pows2[k] + s1*pows1[k] + *p1;
4684 }
4685 }
4686 }
4687 *p2 = *p1;
4688 p1++; p2+=e+1;
4689 }
4690
4691 /* Now go through each column and interpolate the missing rows.
4692 This is done purely within the output grid, since we have already
4693 filled in as much info as we can from the input grid. */
4694
4695 if (exp2==1) return;
4696
4697 for (k=0; k<exp2; k++) {
4698 t=( (float)k ) / ( (float)exp2 ) ;
4699 pows1[k]=t;
4700 pows2[k]=t*t;
4701 pows3[k]=t*pows2[k];
4702 }
4703
4704 p2=g2;
4705
4706 for (jj=0; jj<rows-1; jj++) {
4707 for (ii=0; ii<c; ii++,p2++) {
4708 if (*p2==mval || *(p2+d)==mval) {
4709 for (k=1; k<exp2; k++) *(p2+(c*k))=mval;
4710 } else {
4711 if (jj==0 || *(p2-d)==mval) s1 = *(p2+d) - *p2;
4712 else s1 = ( *(p2+d) - *(p2-d) )*0.5;
4713 if (jj==(rows-2) || *(p2+(2*d))==mval) s2 = *(p2+d) - *p2;
4714 else s2 = ( *(p2+(2*d)) - *p2 )*0.5;
4715 s1*=kurv; s2*=kurv;
4716 a = s1 + 2.0*(*p2) - 2.0*(*(p2+d)) + s2;
4717 b = 3.0*(*(p2+d)) - s2 - 2.0*s1 - 3.0*(*p2);
4718 for (k=1; k<exp2; k++) {
4719 *(p2+(c*k)) = a*pows3[k] + b*pows2[k] + s1*pows1[k] + *p2;
4720 }
4721 }
4722 }
4723 p2+=e;
4724 }
4725 return;
4726 }
4727
4728 /* Rotate a grid. Return rotated grid. */
4729
gaflip(struct gagrid * pgr,struct gacmn * pcm)4730 struct gagrid *gaflip (struct gagrid *pgr, struct gacmn *pcm) {
4731 struct gagrid *pgr2;
4732 float *gr1, *gr2;
4733 int i, j, size;
4734
4735 pgr2 = (struct gagrid *)malloc(sizeof(struct gagrid));
4736 if (pgr2==NULL) {
4737 gaprnt (0,"Memory Allocation Error: gaflip \n");
4738 return (NULL);
4739 }
4740 *pgr2 = *pgr;
4741 size = pgr->isiz * pgr->jsiz;
4742 pgr2->grid = (float *)malloc(size*sizeof(float));
4743 if (pgr2->grid == NULL) {
4744 gaprnt (0,"Memory Allocation Error: gaflip \n");
4745 free (pgr2);
4746 return (NULL);
4747 }
4748
4749 gr1 = pgr->grid;
4750 for (j=0; j<pgr->jsiz; j++) {
4751 gr2 = pgr2->grid + j;
4752 for (i=0; i<pgr->isiz; i++) {
4753 *gr2 = *gr1;
4754 gr1++;
4755 gr2 += pgr->jsiz;
4756 }
4757 }
4758
4759 pgr2->idim = pgr->jdim;
4760 pgr2->jdim = pgr->idim;
4761 pgr2->iwrld = pgr->jwrld; pgr2->jwrld = pgr->iwrld;
4762 pgr2->isiz = pgr->jsiz;
4763 pgr2->jsiz = pgr->isiz;
4764 pgr2->igrab = pgr->jgrab;
4765 pgr2->jgrab = pgr->igrab;
4766 pgr2->ivals = pgr->jvals;
4767 pgr2->jvals = pgr->ivals;
4768 pgr2->ilinr = pgr->jlinr;
4769 pgr2->jlinr = pgr->ilinr;
4770
4771 /* Add new grid to list to be freed later */
4772
4773 i = pcm->relnum;
4774 pcm->type[i] = 1;
4775 pcm->result[i].pgr = pgr2;
4776 pcm->relnum++;
4777
4778 return (pgr2);
4779 }
4780
4781 /* Figure out appropriate date/time increment for time axis. */
4782
gatinc(struct gacmn * pcm,struct dt * tstrt,struct dt * tincr)4783 int gatinc (struct gacmn *pcm, struct dt *tstrt, struct dt *tincr) {
4784 int tdif;
4785 float y1,y2,c1,c2,c3;
4786 struct dt twrk,temp;
4787
4788 tincr->yr = 0L;
4789 tincr->mo = 0L;
4790 tincr->dy = 0L;
4791 tincr->hr = 0L;
4792 tincr->mn = 0L;
4793 twrk.yr = 0L;
4794 twrk.mo = 0L;
4795 twrk.dy = 0L;
4796 twrk.hr = 0L;
4797 twrk.mn = 0L;
4798
4799 /* Check for a time period that covers a period of years. */
4800
4801 if (pcm->tmax.yr-pcm->tmin.yr>4) {
4802 y1 = pcm->tmin.yr;
4803 y2 = pcm->tmax.yr;
4804 c1 = 0.0;
4805 gacsel (y1,y2,&c1,&c2,&c3);
4806 tincr->yr = (int)(c1+0.5);
4807 if (tincr->yr==3) tincr->yr=5;
4808 goto cont;
4809 }
4810
4811 /* Get time difference in minutes */
4812
4813 tdif = timdif (&(pcm->tmin),&(pcm->tmax));
4814
4815 /* Set time increment based on different time differences.
4816 The test is entirely arbitrary. */
4817
4818 if (tdif>1576800L) tincr->mo = 6L;
4819 else if (tdif>788400L) tincr->mo = 3L;
4820 else if (tdif>262800L) tincr->mo = 1L;
4821 else if (tdif>129600L) tincr->dy = 15L;
4822 else if (tdif>42000L) tincr->dy = 5L;
4823 else if (tdif>14399L) tincr->dy = 2L;
4824 else if (tdif>7199L) tincr->dy = 1L;
4825 else if (tdif>3599L) tincr->hr = 12L;
4826 else if (tdif>1799L) tincr->hr = 6L;
4827 else if (tdif>899L) tincr->hr = 3L;
4828 else if (tdif>299L) tincr->hr = 1L;
4829 else if (tdif>149L) tincr->mn = 30L;
4830 else if (tdif>74L) tincr->mn = 15L;
4831 else if (tdif>24L) tincr->mn = 5L;
4832 else if (tdif>9L) tincr->mn = 2L;
4833 else tincr->mn = 1L;
4834
4835 /* Round tmin to get the correct starting time for this increment.
4836 Note that this algorithm assumes that only the above increment
4837 settings are possible. So you can change the range tests
4838 (tdiff>something), but do not change the possible increments. */
4839
4840 cont:
4841 *tstrt = pcm->tmin;
4842
4843 if (tincr->mn>0L) {
4844 tdif = tincr->mn*(tstrt->mn/tincr->mn);
4845 if (tdif!=tstrt->mn) tstrt->mn = tdif+tincr->mn;
4846 if (tstrt->mn>59) {
4847 tstrt->mn = 0L;
4848 temp = twrk;
4849 temp.hr = 1L;
4850 timadd (tstrt,&temp);
4851 *tstrt = temp;
4852 }
4853 return(5);
4854 }
4855 if (tstrt->mn>0L) {
4856 tstrt->mn = 0L;
4857 temp = twrk;
4858 temp.hr = 1L;
4859 timadd (tstrt,&temp);
4860 *tstrt = temp;
4861 }
4862
4863 if (tincr->hr>0L) {
4864 tdif = tincr->hr*(tstrt->hr/tincr->hr);
4865 if (tdif!=tstrt->hr) tstrt->hr = tdif+tincr->hr;
4866 if (tstrt->hr>23) {
4867 tstrt->hr = 0L;
4868 temp = twrk;
4869 temp.dy = 1L;
4870 timadd (tstrt,&temp);
4871 *tstrt = temp;
4872 }
4873 return (4);
4874 }
4875 if (tstrt->hr>0L) {
4876 tstrt->hr = 0L;
4877 temp = twrk;
4878 temp.dy = 1L;
4879 timadd (tstrt,&temp);
4880 *tstrt = temp;
4881 }
4882
4883 if (tincr->dy>0L) {
4884 tdif = 1L+tincr->dy*((tstrt->dy-1L)/tincr->dy);
4885 if (tdif!=tstrt->dy) {
4886 tstrt->dy = tdif+tincr->dy;
4887 if (tstrt->dy==29L || tstrt->dy==31L) {
4888 tstrt->dy = 1L;
4889 temp = twrk;
4890 temp.mo = 1L;
4891 timadd (tstrt,&temp);
4892 *tstrt = temp;
4893 }
4894 }
4895 return (3);
4896 }
4897 if (tstrt->dy>1L) {
4898 tstrt->dy = 1L;
4899 temp = twrk;
4900 temp.mo = 1L;
4901 timadd (tstrt,&temp);
4902 *tstrt = temp;
4903 }
4904
4905 if (tincr->mo>0L) {
4906 tdif = 1L+tincr->mo*((tstrt->mo-1L)/tincr->mo);
4907 if (tdif!=tstrt->mo) tstrt->mo = tdif+tincr->mo;
4908 if (tstrt->mo>12) {
4909 tstrt->mo = 1L;
4910 temp = twrk;
4911 temp.yr = 1L;
4912 timadd (tstrt,&temp);
4913 *tstrt = temp;
4914 }
4915 return (2);
4916 }
4917 if (tstrt->mo>1L) {
4918 tstrt->mo = 1L;
4919 temp = twrk;
4920 temp.yr = 1L;
4921 timadd (tstrt,&temp);
4922 *tstrt = temp;
4923 }
4924
4925 tdif = tincr->yr*(tstrt->yr/tincr->yr);
4926 if (tdif!=tstrt->yr) tstrt->yr = tdif+tincr->yr;
4927 return(1);
4928 }
4929
4930 /* Plot lat/lon lines when polar stereographic plots are done */
4931
gampax(struct gacmn * pcm)4932 void gampax (struct gacmn *pcm) {
4933 float x1,y1,x2,y2,v,s,plincr,lndif,ln,lt,cs;
4934 float lnmin,lnmax,ltmin,ltmax,lnincr,ltincr,lnstrt,lnend,ltstrt,ltend;
4935
4936 if (!pcm->grflag) return;
4937
4938 if (pcm->mproj==5) {
4939 lnmin = pcm->dmin[0]; lnmax = pcm->dmax[0];
4940 ltmin = pcm->dmin[1]; ltmax = pcm->dmax[1];
4941 for (ln=lnmin; ln<lnmax+10.0; ln+=45.0) {
4942 if (ln<lnmin+10.0 || ln>lnmax-10.0) {
4943 gxstyl (1);
4944 gxcolr (1);
4945 gxwide (5);
4946 } else {
4947 gxstyl (pcm->grstyl);
4948 gxcolr (pcm->grcolr);
4949 gxwide (1);
4950 }
4951 gxconv (ln,ltmin,&x1,&y1,2);
4952 gxplot (x1,y1,3);
4953 for (lt=ltmin; lt<ltmax+1.0; lt+=2.0) {
4954 gxconv (ln,lt,&x1,&y1,2);
4955 gxplot (x1,y1,2);
4956 }
4957 }
4958 for (lt=ltmin; lt<ltmax+10.0; lt+=30.0) {
4959 if (lt<ltmin+10.0 || lt>ltmax-10.0) {
4960 gxstyl (1);
4961 gxcolr (1);
4962 gxwide (5);
4963 } else {
4964 gxstyl (pcm->grstyl);
4965 gxcolr (pcm->grcolr);
4966 gxwide (1);
4967 }
4968 gxconv (lnmin,lt,&x1,&y1,2);
4969 gxplot (x1,y1,3);
4970 for (ln=lnmin; ln<lnmax+1.0; ln+=2.0) {
4971 gxconv (ln,lt,&x1,&y1,2);
4972 gxplot (x1,y1,2);
4973 }
4974 }
4975 if (fabs(lnmin+180.0)<1.0 && fabs(lnmax-180.0)<1.0 &&
4976 fabs(ltmin+90.0)<1.0 && fabs(ltmax-90.0)<1.0) {
4977 cs = 0.10;
4978 gxconv (-180.0,90.0,&x1,&y1,2);
4979 gxchpl ("180",3,x1-cs*1.5,y1+cs*0.6,cs*1.1,cs,0.0);
4980 gxconv (-90.0,90.0,&x1,&y1,2);
4981 gxchpl ("90W",3,x1-cs*1.5,y1+cs*0.6,cs*1.1,cs,0.0);
4982 gxconv (0.0,90.0,&x1,&y1,2);
4983 gxchpl ("0",1,x1-cs*0.5,y1+cs*0.6,cs*1.1,cs,0.0);
4984 gxconv (90.0,90.0,&x1,&y1,2);
4985 gxchpl ("90E",3,x1-cs*1.5,y1+cs*0.6,cs*1.1,cs,0.0);
4986 gxconv (180.0,90.0,&x1,&y1,2);
4987 gxchpl ("180",3,x1-cs*1.5,y1+cs*0.6,cs*1.1,cs,0.0);
4988 gxconv (-180.0,-90.0,&x1,&y1,2);
4989 gxchpl ("180",3,x1-cs*1.5,y1-cs*1.6,cs*1.1,cs,0.0);
4990 gxconv (-90.0,-90.0,&x1,&y1,2);
4991 gxchpl ("90W",3,x1-cs*1.5,y1-cs*1.6,cs*1.1,cs,0.0);
4992 gxconv (0.0,-90.0,&x1,&y1,2);
4993 gxchpl ("0",1,x1-cs*0.5,y1-cs*1.6,cs*1.1,cs,0.0);
4994 gxconv (90.0,-90.0,&x1,&y1,2);
4995 gxchpl ("90E",3,x1-cs*1.5,y1-cs*1.6,cs*1.1,cs,0.0);
4996 gxconv (180.0,-90.0,&x1,&y1,2);
4997 gxchpl ("180",3,x1-cs*1.5,y1-cs*1.6,cs*1.1,cs,0.0);
4998
4999 gxconv (-180.0,-60.0,&x1,&y1,2);
5000 gxchpl ("60S",3,x1-cs*4.5,y1-cs*0.55,cs*1.1,cs,0.0);
5001 gxconv (-180.0,-30.0,&x1,&y1,2);
5002 gxchpl ("30S",3,x1-cs*4.0,y1-cs*0.55,cs*1.1,cs,0.0);
5003 gxconv (-180.0,0.0,&x1,&y1,2);
5004 gxchpl ("EQ",2,x1-cs*2.7,y1-cs*0.55,cs*1.1,cs,0.0);
5005 gxconv (-180.0,60.0,&x1,&y1,2);
5006 gxchpl ("60N",3,x1-cs*4.5,y1-cs*0.55,cs*1.1,cs,0.0);
5007 gxconv (-180.0,30.0,&x1,&y1,2);
5008 gxchpl ("30N",3,x1-cs*4.0,y1-cs*0.55,cs*1.1,cs,0.0);
5009
5010 gxconv (180.0,-60.0,&x1,&y1,2);
5011 gxchpl ("60S",3,x1+cs*1.5,y1-cs*0.55,cs*1.1,cs,0.0);
5012 gxconv (180.0,-30.0,&x1,&y1,2);
5013 gxchpl ("30S",3,x1+cs*1.0,y1-cs*0.55,cs*1.1,cs,0.0);
5014 gxconv (180.0,0.0,&x1,&y1,2);
5015 gxchpl ("EQ",2,x1+cs*0.7,y1-cs*0.55,cs*1.1,cs,0.0);
5016 gxconv (180.0,60.0,&x1,&y1,2);
5017 gxchpl ("60N",3,x1+cs*1.5,y1-cs*0.55,cs*1.1,cs,0.0);
5018 gxconv (180.0,30.0,&x1,&y1,2);
5019 gxchpl ("30N",3,x1+cs*1.0,y1-cs*0.55,cs*1.1,cs,0.0);
5020
5021 }
5022 return;
5023 } /*mf eeeeeeeeee end of robinson projection case */
5024
5025 /* Choose grid interval for longitude */
5026
5027 lnincr=0.0;
5028 lnmin = pcm->dmin[0];
5029 lnmax = pcm->dmax[0];
5030 gacsel (lnmin,lnmax,&lnincr,&lnstrt,&lnend);
5031 /* printf("qqq %f %f %f %f %f\n",lnmin,lnmax,lnincr,lnstrt,lnend); */
5032 if (lnincr==0.0) {
5033 gaprnt (0,"gaaxis internal logic check 25\n");
5034 return;
5035 }
5036 if (lnincr>24.9) lnincr=30.0;
5037 else if (lnincr>14.5) lnincr=20.0;
5038 else if (lnincr>10.0) lnincr=10.0;
5039 gacsel (lnmin,lnmax,&lnincr,&lnstrt,&lnend);
5040 lndif = lnmax - lnmin;
5041
5042 if(pcm->xlint!=0.0) lnincr=pcm->xlint; /*mf 960402 set the lon increment from xlint */
5043
5044 /* Choose grid interval for latitude */
5045
5046 ltincr=0.0;
5047 ltmin = pcm->dmin[1];
5048 ltmax = pcm->dmax[1];
5049 gacsel (ltmin,ltmax,<incr,<strt,<end);
5050 if (ltincr==0.0) {
5051 gaprnt (0,"gaaxis internal logic check 25\n");
5052 return;
5053 }
5054 if (lndif>300.0) {
5055 if (ltincr>9.9) ltincr = 20.0;
5056 else if (ltincr>4.9) ltincr = 10.0;
5057 else if (ltincr>1.9) ltincr = 5.0;
5058 else if (ltincr>0.8) ltincr = 2.0;
5059 else if (ltincr>0.4) ltincr = 1.0;
5060 else if (ltincr>0.2) ltincr = 0.5;
5061 } else {
5062 if (ltincr>14.9) ltincr = 20.0;
5063 else if (ltincr>7.9) ltincr = 10.0;
5064 else if (ltincr>2.9) ltincr = 5.0;
5065 else if (ltincr>1.4) ltincr = 2.0;
5066 else if (ltincr>0.6) ltincr = 1.0;
5067 else if (ltincr>0.2) ltincr = 0.5;
5068 }
5069
5070 if(pcm->ylint!=0.0) ltincr=pcm->ylint; /*mf 960402 set the lat increment from ylint */
5071
5072 if (ltstrt<-89.9) {
5073 ltstrt = ltstrt + ltincr;
5074 ltmin = ltstrt;
5075 }
5076 if (ltend>89.9) {
5077 ltend = ltend - ltincr;
5078 ltmax = ltend;
5079 }
5080
5081 gxstyl (pcm->grstyl);
5082 gxcolr (pcm->grcolr);
5083 gxwide (1);
5084 gxclip (pcm->xsiz1,pcm->xsiz2,pcm->ysiz1,pcm->ysiz2);
5085 for (v=lnstrt; v<lnend+lnincr*0.5; v+=lnincr) {
5086 gxconv (v,ltmin,&x1,&y1,2);
5087 gxconv (v,ltmax,&x2,&y2,2);
5088 gxplot (x1,y1,3);
5089 gxplot (x2,y2,2);
5090 }
5091
5092 if (lndif>180.0) plincr = lndif/100.0;
5093 else if (lndif>60.0) plincr = lndif/50.0;
5094 else plincr = lndif/25.0;
5095 for (v=ltstrt; v<ltend+ltincr*0.5; v+=ltincr) {
5096 gxconv (lnmin,v,&x1,&y1,2);
5097 gxplot (x1,y1,3);
5098 for (s=lnmin+plincr; s<lnmax+plincr*0.5; s+=plincr) {
5099 gxconv (s,v,&x1,&y1,2);
5100 gxplot (x1,y1,2);
5101 }
5102 }
5103
5104 /* Plot circular frame if user requested such a thing */
5105
5106 if (pcm->frame==2) {
5107 gxcolr(pcm->anncol);
5108 gxwide (pcm->annthk-3);
5109 gxstyl (1);
5110 if (pcm->mproj==3) v = ltmin;
5111 else v = ltmax;
5112 gxconv (lnmin,v,&x1,&y1,2);
5113 gxplot (x1,y1,3);
5114 for (s=lnmin+plincr; s<lnmax+plincr*0.5; s+=plincr) {
5115 gxconv (s,v,&x1,&y1,2);
5116 gxplot (x1,y1,2);
5117 }
5118 }
5119
5120 gxclip (0.0, pcm->xsiz, 0.0, pcm->ysiz);
5121 gxstyl(1);
5122 }
5123
5124 /* Plot weather symbol */
5125
wxsym(int type,float xpos,float ypos,float scale,int colr,int * wxcols)5126 void wxsym (int type, float xpos, float ypos, float scale, int colr,
5127 int *wxcols) {
5128 int ioff,len,i,part,icol;
5129 float x,y;
5130
5131 wxymin = 9e30;
5132 wxymax = -9e30;
5133 if (type<1||type>43) return;
5134 ioff = strt[type-1]-1;
5135 len = scnt[type-1];
5136 for (i=0; i<len; i++) {
5137 part = stype[i+ioff]-1;
5138 if (colr<0) {
5139 if (type==25) icol=0;
5140 else if (type==14 || type==15) icol=1;
5141 else if (type==19 || type==20) icol=2;
5142 else if (type==31 || type==32) icol=0;
5143 else icol = scol[part];
5144 gxcolr (*(wxcols+icol));
5145 } else gxcolr (colr);
5146 x = xpos + sxpos[i+ioff]*scale;
5147 y = ypos + sypos[i+ioff]*scale;
5148 wxprim (part, x, y, scale);
5149 }
5150 if (type==40 || type==42) gxmark (2,x,y,scale*0.4);
5151 if (type==41 || type==43) gxmark (3,x,y,scale*0.4);
5152 }
5153
5154
5155 /* Plot primitive wx symbol */
5156
wxprim(int part,float xpos,float ypos,float scale)5157 void wxprim (int part, float xpos, float ypos, float scale) {
5158 int i,len,pos;
5159 float x,y;
5160 float xy[40];
5161
5162 /* xxxx */
5163 len = slen[part];
5164 pos = soff[part]-1;
5165 for (i=0; i<len; i++) {
5166 x = xpts[i+pos];
5167 y = ypts[i+pos];
5168 x = x*scale + xpos;
5169 y = y*scale + ypos;
5170 if (part==1 | part==2) {
5171 xy[i*2] = x;
5172 xy[i*2+1] = y;
5173 } else {
5174 gxplot (x,y,spens[i+pos]);
5175 }
5176 if (y<wxymin) wxymin = y;
5177 if (y>wxymax) wxymax = y;
5178 }
5179 if (part==1 | part==2) gxfill(xy,len);
5180 }
5181
gagsav(int type,struct gacmn * pcm,struct gagrid * pgr)5182 void gagsav (int type, struct gacmn *pcm, struct gagrid *pgr) {
5183 pcm->lastgx = type;
5184 }
5185
5186 /* Log axis scaling */
5187
galnx(float xin,float yin,float * xout,float * yout)5188 void galnx (float xin, float yin, float *xout, float *yout) {
5189 *xout = log(xin);
5190 *yout = yin;
5191 }
5192
galny(float xin,float yin,float * xout,float * yout)5193 void galny (float xin, float yin, float *xout, float *yout) {
5194 *xout = xin;
5195 *yout = log(yin);
5196 }
5197
gaalnx(float xin,float yin,float * xout,float * yout)5198 void gaalnx (float xin, float yin, float *xout, float *yout) {
5199 *xout = pow(2.71829,xin);
5200 *yout = yin;
5201 }
5202
gaalny(float xin,float yin,float * xout,float * yout)5203 void gaalny (float xin, float yin, float *xout, float *yout) {
5204 *xout = xin;
5205 *yout = pow(2.71829,yin);
5206 }
5207
5208 /* cos lat scaling */
5209
gaclx(float xin,float yin,float * xout,float * yout)5210 void gaclx (float xin, float yin, float *xout, float *yout) {
5211 *xout = sin(xin*3.1416/180.0);
5212 *yout = yin;
5213 }
5214
gacly(float xin,float yin,float * xout,float * yout)5215 void gacly (float xin, float yin, float *xout, float *yout) {
5216 *xout = xin;
5217 *yout = sin(yin*3.1416/180.0);
5218 }
5219
gaaclx(float xin,float yin,float * xout,float * yout)5220 void gaaclx (float xin, float yin, float *xout, float *yout) {
5221 *xout = asin(xin)*180.0/3.1416;
5222 *yout = yin;
5223 }
5224
gaacly(float xin,float yin,float * xout,float * yout)5225 void gaacly (float xin, float yin, float *xout, float *yout) {
5226 *xout = xin;
5227 *yout = asin(yin)*180.0/3.1416;
5228 }
5229
5230
5231 /* Grid fill -- using color ranges as in shaded contours */
5232
gagfil(float * r,int is,int js,float * vs,int * clrs,int lvs,float u)5233 void gagfil ( float *r, int is, int js, float *vs, int *clrs, int lvs,
5234 float u) {
5235 int i,j,k,flag;
5236 float x,y,xhi,yhi;
5237 float *v,xybox[10];
5238
5239 v = r;
5240 for (j=1; j<=js; j++) {
5241 for (i=1; i<=is; i++) {
5242 if (*v != u) {
5243 flag = 1;
5244 for (k=0; k<lvs-1; k++) {
5245 if (*v>*(vs+k) && *v<=*(vs+k+1)) {
5246 gxcolr(*(clrs+k));
5247 gxconv ((float)(i)-0.5,(float)(j)-0.5,&x,&y,3);
5248 xybox[0] = x; xybox[1] = y;
5249 gxconv ((float)(i)+0.5,(float)(j)-0.5,&x,&y,3);
5250 xybox[2] = x; xybox[3] = y;
5251 gxconv ((float)(i)+0.5,(float)(j)+0.5,&x,&y,3);
5252 xybox[4] = x; xybox[5] = y;
5253 gxconv ((float)(i)-0.5,(float)(j)+0.5,&x,&y,3);
5254 xybox[6] = x; xybox[7] = y;
5255 gxconv ((float)(i)-0.5,(float)(j)-0.5,&x,&y,3);
5256 xybox[8] = x; xybox[9] = y;
5257 gxfill (xybox,5);
5258 flag = 0;
5259 break;
5260 }
5261 }
5262 if (flag) {
5263 gxcolr(*(clrs+lvs-1));
5264 gxconv ((float)(i)-0.5,(float)(j)-0.5,&x,&y,3);
5265 xybox[0] = x; xybox[1] = y;
5266 gxconv ((float)(i)+0.5,(float)(j)-0.5,&x,&y,3);
5267 xybox[2] = x; xybox[3] = y;
5268 gxconv ((float)(i)+0.5,(float)(j)+0.5,&x,&y,3);
5269 xybox[4] = x; xybox[5] = y;
5270 gxconv ((float)(i)-0.5,(float)(j)+0.5,&x,&y,3);
5271 xybox[6] = x; xybox[7] = y;
5272 gxconv ((float)(i)-0.5,(float)(j)-0.5,&x,&y,3);
5273 xybox[8] = x; xybox[9] = y;
5274 gxfill (xybox,5);
5275 }
5276 }
5277 v++;
5278 } }
5279 }
5280
gafram(struct gacmn * pcm)5281 void gafram(struct gacmn *pcm) {
5282
5283 if (pcm->frame==0) return;
5284 if (pcm->frame==2 && pcm->mproj>2) return;
5285 if (pcm->mproj>4) return;
5286
5287 gxcolr (pcm->anncol);
5288 gxwide (pcm->annthk);
5289 gxstyl (1);
5290 gxplot (pcm->xsiz1,pcm->ysiz1,3);
5291 gxplot (pcm->xsiz2,pcm->ysiz1,2);
5292 gxplot (pcm->xsiz2,pcm->ysiz2,2);
5293 gxplot (pcm->xsiz1,pcm->ysiz2,2);
5294 gxplot (pcm->xsiz1,pcm->ysiz1,2);
5295 }
5296
gaaxpl(struct gacmn * pcm,int idim,int jdim)5297 void gaaxpl (struct gacmn *pcm, int idim, int jdim) {
5298
5299 if (idim==0 && jdim==1 && pcm->mproj>2) gampax(pcm);
5300 else {
5301 gaaxis(1,pcm,idim);
5302 gaaxis(0,pcm,jdim);
5303 }
5304 }
5305
5306 /* Select colors and intervals for colorized items such as
5307 colorized vectors, streamlines, barbs, scatter, etc. */
5308
gaselc(struct gacmn * pcm,float rmin,float rmax)5309 void gaselc (struct gacmn *pcm, float rmin, float rmax) {
5310 float rdif,w1,norml,cint,cmin,cmax,t1,t2,fact,rr,rrb,smin,ci,ccnt;
5311 int i,ii,irb,power,lcnt;
5312
5313 rdif = rmax - rmin;
5314 if (rdif==0.0) {
5315 pcm->shdcnt = 0;
5316 return;
5317 }
5318
5319 if (pcm->rbflg) irb = pcm->rbflg - 1;
5320 else irb = 12;
5321 rrb = irb+1;
5322
5323 /* Choose a contour interval if one is needed. We choose
5324 an interval that maximizes the number of colors displayed,
5325 rather than an interval that maximizes the attribute of
5326 being a 'nice' interval */
5327
5328 if (!pcm->cflag) {
5329 if (pcm->cint==0.0) {
5330 if (pcm->rbflg) rdif = rdif/((float)(pcm->rbflg)-0.5);
5331 else rdif = rdif/12.5;
5332 w1 = floor(log10(rdif));
5333 w1 = pow(10.0,w1);
5334
5335 norml = rdif/w1; /* normalized contour interval */
5336 if (norml<2.0) fact = 10.0;
5337 else if (norml>=2.0 && norml<4.0) fact = 5.0;
5338 else if (norml>=4.0 && norml<7.0) fact = 2.0;
5339 else fact = 1.0;
5340 norml*=fact;
5341 i = (int)(norml+0.5);
5342 cint = (float)i * w1 / fact;
5343 } else cint = pcm->cint;
5344
5345 cmin = cint * ceil(rmin/cint);
5346 cmax = cint * floor(rmax/cint);
5347
5348 /* Check for interval being below machine epsilon for these
5349 values */
5350
5351 t1 = rmin + cint;
5352 t2 = rmax - cint;
5353 if (t1 == rmin || t2 == rmax) {
5354 pcm->shdcnt = 0;
5355 return;
5356 }
5357
5358 /* Figure out actual contour levels and colors, based on
5359 possible user settings such as cmin, etc. */
5360
5361 if (pcm->rainmn==0.0 && pcm->rainmx==0.0) {
5362 pcm->rainmn = cmin;
5363 pcm->rainmx = cmax;
5364 }
5365 lcnt=0;
5366 smin = pcm->rainmn - cint;
5367 for (rr=cmin-cint;rr<=cmax+(cint/2.0);rr+=cint) {
5368 if (rr<0.0) i = (int)((rr/cint)-0.1);
5369 else i = (int)((rr/cint)+0.1);
5370 rr = (float)i * cint;
5371 pcm->shdlvs[lcnt] = rr;
5372 pcm->shdcls[lcnt] = 1;
5373 if (rr<pcm->cmin) pcm->shdcls[lcnt] = -1;
5374 if (rr+cint>pcm->cmax) pcm->shdcls[lcnt] = -1;
5375 if (pcm->blkflg && rr<pcm->blkmax && rr+cint>pcm->blkmin)
5376 pcm->shdcls[lcnt] = -1;
5377 if (pcm->shdcls[lcnt]==1) {
5378 ii = (int)(rrb*(rr-smin)/(pcm->rainmx-smin));
5379 if (ii>irb) ii=irb;
5380 if (ii<0) ii=0;
5381 if (pcm->rbflg) pcm->shdcls[lcnt] = pcm->rbcols[ii];
5382 else pcm->shdcls[lcnt] = rcols[ii];
5383 }
5384 if (lcnt<1 || pcm->shdcls[lcnt]!=pcm->shdcls[lcnt-1]) lcnt++;
5385 }
5386 pcm->shdlvs[0] -= cint*10.0;
5387 pcm->shdcnt = lcnt;
5388
5389 /* User has specified actual contour levels. Just use those */
5390
5391 } else {
5392
5393 if (rmin<pcm->clevs[0]) pcm->shdlvs[0] = rmin-100.0;
5394 else pcm->shdlvs[0] = pcm->clevs[0]-100.0;
5395 ccnt = rrb/(float)(pcm->cflag+1);
5396 ci = ccnt/2.0;
5397 ii = (int)ci;
5398 if (pcm->ccflg>0) pcm->shdcls[0] = pcm->ccols[0];
5399 else {
5400 if (pcm->rbflg>0) pcm->shdcls[0] = pcm->rbcols[ii];
5401 else pcm->shdcls[0] = rcols[ii];
5402 }
5403 lcnt = 1;
5404 for (i=0; i<pcm->cflag; i++) {
5405 pcm->shdlvs[lcnt] = pcm->clevs[i];
5406 if (pcm->ccflg>0) {
5407 if (i+1>=pcm->ccflg) pcm->shdcls[lcnt]=15;
5408 else pcm->shdcls[lcnt]=pcm->ccols[i+1];
5409 } else {
5410 ci += ccnt;
5411 ii = (int)ci;
5412 if (ii>irb) ii=irb;
5413 if (ii<0) ii=0;
5414 if (pcm->rbflg>0) pcm->shdcls[lcnt] = pcm->rbcols[ii];
5415 else pcm->shdcls[lcnt] = rcols[ii];
5416 }
5417 if (lcnt<1 || pcm->shdcls[lcnt]!=pcm->shdcls[lcnt-1]) lcnt++;
5418 }
5419 pcm->shdcnt = lcnt;
5420 }
5421 }
5422
5423 /* Given a shade value, return the relevent color */
5424
gashdc(struct gacmn * pcm,float val)5425 int gashdc (struct gacmn *pcm, float val) {
5426 int i;
5427
5428 if (pcm->shdcnt==0) return(1);
5429 if (pcm->shdcnt==1) return(pcm->shdcls[0]);
5430 if (val<pcm->shdlvs[1]) return(pcm->shdcls[0]);
5431 for (i=1; i<pcm->shdcnt-1; i++) {
5432 if (val>=pcm->shdlvs[i] && val<pcm->shdlvs[i+1])
5433 return(pcm->shdcls[i]);
5434 }
5435 return(pcm->shdcls[pcm->shdcnt-1]);
5436 }
5437
5438
5439
5440
5441 void gatmlb (struct gacmn *); /*mf --- local protype --- mf*/
5442
5443
5444 struct cxclock {
5445 time_t year ;
5446 time_t month ;
5447 time_t date ;
5448 time_t hour ;
5449 time_t minute ;
5450 time_t second ;
5451 char timezone[32] ;
5452 char clockenv[64] ;
5453 int julian_day ;
5454 time_t epoch_time_in_sec ;
5455 } ;
5456
5457 /* global variable of the time object */
5458 static struct cxclock timeobj;
5459
5460
get_tm_struct(time_t t)5461 void get_tm_struct(time_t t)
5462 {
5463 struct tm *ptr_tm ;
5464 char buf[10] ;
5465 time_t sec ;
5466 sec=t ;
5467
5468 /* get tm struct from seconds */
5469 ptr_tm=localtime(&sec) ;
5470
5471 if ( ! ptr_tm ) {
5472 fprintf(stderr,
5473 "error in gmtime function,errno = %d\n",errno) ;
5474 exit(-1) ;
5475 }
5476
5477 /* copy tm struct to cxclock private members */
5478
5479 timeobj.epoch_time_in_sec=sec ;
5480 timeobj.year=ptr_tm->tm_year + 1900 ;
5481 timeobj.month=(ptr_tm->tm_mon) + 1 ;
5482 timeobj.date=ptr_tm->tm_mday ;
5483 timeobj.hour=ptr_tm->tm_hour ;
5484 timeobj.minute=ptr_tm->tm_min ;
5485 timeobj.second=ptr_tm->tm_sec ;
5486
5487 /* get julian date */
5488 strftime(buf, sizeof(buf),"%j", ptr_tm) ;
5489 timeobj.julian_day=atoi(buf) ;
5490
5491 }
5492
sys_time()5493 void sys_time()
5494 {
5495 time_t sec ;
5496 time(&sec) ; /* get the GMT offset from the local machine */
5497 if (sec == -1 ) {
5498 fprintf(stderr,
5499 "error in time function,errno = %d\n"
5500 ,errno) ;
5501 exit(-1) ;
5502 }
5503
5504 /* Set tm struct using new epoch time */
5505 get_tm_struct(sec) ;
5506 return;
5507 }
5508
5509
gatmlb(struct gacmn * pcm)5510 void gatmlb (struct gacmn *pcm) {
5511
5512 int i,vcnt;
5513 char dtgstr[32] ;
5514
5515 vcnt = 0;
5516 for (i=0; i<4; i++) if (pcm->vdim[i]) vcnt++;
5517
5518 /* -- only plot if 2-D or less, i.e., turn off on looping */
5519
5520 if (vcnt<=2) {
5521
5522 sys_time();
5523
5524 sprintf(dtgstr,"%04d-%02d-%02d-%02d:%02d\0",
5525 timeobj.year,timeobj.month,
5526 timeobj.date,timeobj.hour,timeobj.minute) ;
5527
5528 gxchpl(dtgstr,strlen(dtgstr),pcm->pxsize-1.36,0.05,0.1,0.08,0.0);
5529
5530 }
5531
5532 }
5533
5534