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 /* Originally authored by B. Doty */
7 /* kk 020624 --- change for 64bit seek K.Komine */
8 #include <stdio.h>
9 #include <math.h>
10 #include <limits.h>
11
12 /*
13 * Include ./configure's header file
14 */
15 #ifdef HAVE_CONFIG_H
16
17 #include <config.h>
18 #if READLINE == 1
19 #include <sys/types.h>
20 #include <sys/file.h>
21 #include <sys/stat.h>
22 #include <sys/errno.h>
23 #include "readline/readline.h"
24 #include "readline/history.h"
25 #endif
26
27 /* If autoconfed, only include malloc.h when it's presen */
28 #ifdef HAVE_MALLOC_H
29 #include <malloc.h>
30 #endif
31
32 #else /* undef HAVE_CONFIG_H */
33 #if READLINE == 1
34 #include <sys/types.h>
35 #include <sys/file.h>
36 #include <sys/stat.h>
37 #include <sys/errno.h>
38 #include "readline.h"
39 #include "history.h"
40 #endif
41
42 #include <malloc.h>
43
44 #endif /* HAVE_CONFIG_H */
45
46 #include "grads.h"
47 #include "gx.h"
48
49 /*mf 961205 --- expose Mike Fiorino's global struct to these routines for 365 day calandars mf*/
50 /*-------- declare here vice grads.c because this is a common routine with LATS ------*/
51 struct gamfcmn mfcmn;
52 /*mf 961205 --- expose Mike Fiorino's global struct to these routines for 365 day calandars mf*/
53
54 static char pout[256]; /* Build Error msgs here */
55
56 /*ams 971228
57 On MS-DOS with XLIBEMU we handle user input differently.
58 See pcx11e.c for a nxtcmd() replacement. However, the utilities
59 don't need this new user interface - macro STNDALN takes care
60 of that.
61 ams*/
62
63 #if !defined(XLIBEMU) || defined(STNDALN)
64
65 /* Retrieves the next command from the user. Leading blanks
66 are stripped. The number of characters entered before the
67 CR is returned. */
68
nxtcmd(char * cmd,char * prompt)69 int nxtcmd (char *cmd, char *prompt) {
70 int i,past,cnt;
71 static char *ch;
72
73 printf ("%s ",prompt);
74 past = 0;
75 cnt = 0;
76 while (1) {
77 *cmd = getchar();
78 if (*cmd == EOF) return (-1);
79 if (*cmd == '\n') {
80 cmd++;
81 *cmd = '\0';
82 return (cnt);
83 }
84 if (past || *cmd != ' ') {
85 cmd++; cnt++; past = 1;
86 }
87 }
88 }
89
90 #endif /* XLIBEMU */
91
92 /* Retrieves the next command from the user. Leading blanks
93 are stripped. The number of characters entered before the
94 CR is returned. */
95
96 #if READLINE == 1
nxrdln(char * cmd,char * prompt)97 int nxrdln (char *cmd, char *prompt) {
98 char *ch, *ch2;
99
100 if((ch=readline(prompt)) == NULL) {
101 return(-1);
102 } else {
103 ch2 = ch;
104 while (*ch == ' ') ch++; /* Skip leading blanks */
105 strcpy(cmd, ch);
106 if (*ch) add_history(ch); /* Skip blank lines */
107 free(ch2);
108 }
109 return(strlen(cmd)+1);
110
111 }
112 #endif
113
114 /* Date/Time manipulation routines. Note that these routines
115 are not particularly efficient, thus Date/Time conversions
116 should be kept to a minimum. */
117
118 static int mosiz[13] = {0,31,28,31,30,31,30,31,31,30,31,30,31};
119 static int momn[13] = {0,44640,40320,44640,43200,44640,43200,
120 44640,44640,43200,44640,43200,44640};
121 static int mnacum[13] = {0,0,44640,84960,129600,172800,217440,
122 260640,305280,349920,393120,437760,480960};
123 static int mnacul[13] = {0,0,44640,86400,131040,174240,218880,
124 262080,306720,351360,394560,439200,482400};
125
126 /* Add an offset to a time. Output to dto. */
127
timadd(struct dt * dtim,struct dt * dto)128 void timadd (struct dt *dtim, struct dt *dto) {
129 int i;
130 int cont;
131
132 /* First add months and years. Normalize as needed. */
133
134 dto->mo += dtim->mo;
135 dto->yr += dtim->yr;
136
137 while (dto->mo>12) {
138 dto->mo -= 12;
139 dto->yr++;
140 }
141
142 /* Add minutes, hours, and days directly. Then normalize
143 to days, then normalize extra days to months/years. */
144
145 dto->mn += dtim->mn;
146 dto->hr += dtim->hr;
147 dto->dy += dtim->dy;
148
149 if (dto->mn > 59) {
150 i = dto->mn / 60;
151 dto->hr += i;
152 dto->mn = dto->mn - (i*60);
153 }
154 if (dto->hr > 23) {
155 i = dto->hr / 24;
156 dto->dy += i;
157 dto->hr = dto->hr - (i*24);
158 }
159
160 cont = 1;
161 while (dto->dy > mosiz[dto->mo] && cont) {
162 if (dto->mo==2 && qleap(dto->yr)) {
163 if (dto->dy == 29) cont=0;
164 else {
165 dto->dy -= 29;
166 dto->mo++;
167 }
168 } else {
169 dto->dy -= mosiz[dto->mo];
170 dto->mo++;
171 }
172 while (dto->mo > 12) {dto->mo-=12; dto->yr++;}
173 }
174 }
175
176 /* Subtract an offset from a time. Subtract minutes/hours/days
177 first so that we will exactly reverse the operation of timadd */
178
timsub(struct dt * dtim,struct dt * dto)179 void timsub (struct dt *dtim, struct dt *dto) {
180 int s1,s2;
181
182 /* Subtract minutes, hour, and days directly. Then normalize
183 to days, then normalize deficient days from months/years. */
184
185 dto->mn = dtim->mn - dto->mn;
186 dto->hr = dtim->hr - dto->hr;
187 dto->dy = dtim->dy - dto->dy;
188 s1 = dto->mo; s2 = dto->yr;
189 dto->mo = dtim->mo;
190 dto->yr = dtim->yr;
191
192 while (dto->mn < 0) {dto->mn+=60; dto->hr--;}
193 while (dto->hr < 0) {dto->hr+=24; dto->dy--;}
194
195 while (dto->dy < 1) {
196 dto->mo--;
197 if (dto->mo < 1) {dto->mo=12; dto->yr--;}
198 if (dto->mo==2 && qleap(dto->yr)) dto->dy += 29;
199 else dto->dy += mosiz[dto->mo];
200 }
201
202 /* Now subtract months and years. Normalize as needed. */
203
204 dto->mo = dto->mo - s1;
205 dto->yr = dto->yr - s2;
206
207 while (dto->mo < 1) {dto->mo+=12; dto->yr--;}
208
209 /* Adjust for leaps */
210
211 if (dto->mo==2 && dto->dy==29 && !qleap(dto->yr)) {
212 dto->mo=3; dto->dy=1;
213 }
214 }
215
216 /* Convert from Absolute time (year/month/day/etc.) to grid
217 coordinate. */
218
t2gr(float * vals,struct dt * dtim)219 float t2gr (float *vals, struct dt *dtim) {
220 struct dt stim;
221 int eyear,mins;
222 float val,*moincr,*mnincr;
223 float rdiff;
224
225 /* Get constants associated with this conversion */
226
227 stim.yr = (int)(*vals+0.1);
228 stim.mo = (int)(*(vals+1)+0.1);
229 stim.dy = (int)(*(vals+2)+0.1);
230 stim.hr = (int)(*(vals+3)+0.1);
231 stim.mn = (int)(*(vals+4)+0.1);
232
233 moincr = vals+5;
234 mnincr = vals+6;
235
236 /* If the increment for this conversion is days, hours, or minutes,
237 then we do our calculations in minutes. If the increment is
238 months or years, we do our calculations in months. */
239
240 if (*mnincr>0.1) {
241 mins = timdif(&stim,dtim);
242 rdiff = (float)mins;
243 val = rdiff/(*mnincr);
244 val += 1.0;
245 return (val);
246 } else {
247 eyear = stim.yr;
248 if (stim.yr > dtim->yr) eyear = dtim->yr;
249 rdiff = (((dtim->yr - eyear)*12) + dtim->mo) -
250 (((stim.yr - eyear)*12) + stim.mo);
251 stim.yr = dtim->yr;
252 stim.mo = dtim->mo;
253 mins = timdif(&stim,dtim);
254 if (mins>0) {
255 if (dtim->mo==2 && qleap(dtim->yr) ) {
256 rdiff = rdiff + (((float)mins)/41760.0);
257 } else {
258 rdiff = rdiff + (((float)mins)/((float)momn[dtim->mo]));
259 }
260 }
261 val = rdiff/(*moincr);
262 val += 1.0;
263 return (val);
264 }
265 }
266
267 /* Convert from a t grid coordinate to an absolute time. */
268
gr2t(float * vals,float gr,struct dt * dtim)269 void gr2t (float *vals, float gr, struct dt *dtim) {
270 struct dt stim;
271 float *moincr,*mnincr;
272 float v;
273
274 /* Get constants associated with this conversion */
275 stim.yr = (int)(*vals+0.1);
276 stim.mo = (int)(*(vals+1)+0.1);
277 stim.dy = (int)(*(vals+2)+0.1);
278 stim.hr = (int)(*(vals+3)+0.1);
279 stim.mn = (int)(*(vals+4)+0.1);
280 moincr = vals+5;
281 mnincr = vals+6;
282
283 /* Initialize output time */
284 dtim->yr = 0;
285 dtim->mo = 0;
286 dtim->dy = 0;
287 dtim->hr = 0;
288 dtim->mn = 0;
289
290 /* Do conversion if increment is in minutes. */
291 if (*mnincr>0.1) {
292 v = *mnincr * (gr-1.0);
293 if (v>0.0) v = v + 0.5; /* round */
294 else v = v - 0.5;
295 dtim->mn = (int)v;
296 if (dtim->mn<0) {
297 dtim->mn = -1 * dtim->mn;
298 timsub (&stim,dtim);
299 } else {
300 timadd (&stim,dtim);
301 }
302 return;
303
304 /* Do conversion if increment is in months. Same as for minutes,
305 except special handling is required for partial months. */
306
307 } else {
308 v = *moincr * (gr-1.0);
309 if (v<0.0) dtim->mo = (int)(v-0.9999); /* round (sort of) */
310 else dtim->mo = (int)(v+0.0001);
311 v = v - (float)dtim->mo; /* Get fractional month */
312 if (dtim->mo<0) {
313 dtim->mo = -1 * dtim->mo;
314 timsub (&stim,dtim);
315 } else timadd (&stim,dtim);
316 if (v<0.0001) return; /* if fraction small, return */
317
318 if (dtim->mo==2 && qleap(dtim->yr) ) {
319 v = v * 41760.0;
320 } else {
321 v = v * (float)momn[dtim->mo];
322 }
323 stim = *dtim;
324 dtim->yr = 0;
325 dtim->mo = 0;
326 dtim->dy = 0;
327 dtim->hr = 0;
328 dtim->mn = (int)(v+0.5);
329 timadd (&stim,dtim);
330 return;
331 }
332 }
333
334 /* Calculate the difference between two times and return the
335 difference in minutes. The calculation is time2 - time1, so
336 if time2 is earlier than time1, the result is negative. */
337
timdif(struct dt * dtim1,struct dt * dtim2)338 int timdif (struct dt *dtim1, struct dt *dtim2) {
339 int min1,min2,yr;
340 struct dt *temp;
341 int swap,mo1,mo2;
342
343 swap = 0;
344 if (dtim1->yr > dtim2->yr) {
345 temp = dtim1;
346 dtim1 = dtim2;
347 dtim2 = temp;
348 swap = 1;
349 }
350
351 min1 = 0;
352 min2 = 0;
353
354 yr = dtim1->yr;
355 while (yr < dtim2->yr) {
356 if (qleap(yr)) min2 += 527040L;
357 else min2 += 525600L;
358 yr++;
359 }
360
361 mo1 = dtim1->mo;
362 mo2 = dtim2->mo;
363 if (qleap(dtim1->yr)) {
364 min1 = min1+mnacul[mo1]+(dtim1->dy*1440L)+(dtim1->hr*60L)+dtim1->mn;
365 } else {
366 min1 = min1+mnacum[mo1]+(dtim1->dy*1440L)+(dtim1->hr*60L)+dtim1->mn;
367 }
368 if (qleap(dtim2->yr)) {
369 min2 = min2+mnacul[mo2]+(dtim2->dy*1440L)+(dtim2->hr*60L)+dtim2->mn;
370 } else {
371 min2 = min2+mnacum[mo2]+(dtim2->dy*1440L)+(dtim2->hr*60L)+dtim2->mn;
372 }
373 if (swap) return (min1-min2);
374 else return (min2-min1);
375 }
376
377 /* Test for leap year. Rules are:
378
379 Divisible by 4, it is a leap year, unless....
380 Divisible by 100, it is not a leap year, unless...
381 Divisible by 400, it is a leap year. */
382
qleap(int year)383 int qleap (int year) {
384 int i,y;
385
386 /*mf - disable if 365 day calendar mf*/
387
388 if(mfcmn.cal365 == 1) return(0);
389
390 y = year;
391
392 i = y / 4;
393 i = (i*4) - y;
394 if (i!=0) return (0);
395
396 i = y / 100;
397 i = (i*100) - y;
398 if (i!=0) return (1);
399
400 i = y / 400;
401 i = (i*400) - y;
402 if (i!=0) return (0);
403
404 return (1);
405
406
407 }
408
409 static char *mons[12] = {"jan","feb","mar","apr","may","jun","jul","aug",
410 "sep","oct","nov","dec"};
411
412 /* Parse an absolute date/time value. Format is:
413
414 12:00z 1jan 1989 (jan,feb,mar,apr,may,jun,jul,aug,sep,oct,nov,dec)
415
416 Must have Z or Month abbrev, or value is invalid. 'def' contains
417 higher order missing values (usually from tmin in pst). Lower order
418 values are defaulted to be: dy = 1, hr = 0, mn = 0. */
419
adtprs(char * ch,struct dt * def,struct dt * dtim)420 char *adtprs (char *ch, struct dt *def, struct dt *dtim) {
421 int val,flag,i;
422 char *pos;
423 char monam[5];
424
425 pos = ch;
426
427 dtim->mn = 0;
428 dtim->hr = 0;
429 dtim->dy = 1;
430
431 if (*ch>='0' && *ch<='9') {
432 flag = 0;
433 ch = intprs (ch,&val);
434 if (*ch == ':' || tolower(*ch) == 'z') {
435 if (val>23) {
436 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
437 sprintf (pout," Hour = %i -- greater than 23\n",val);
438 gaprnt (0,pout);
439 return (NULL);
440 }
441 dtim->hr = val;
442 if (*ch == ':') {
443 ch++;
444 if (*ch>='0' && *ch<='9') {
445 ch = intprs (ch,&val);
446 if (val>59) {
447 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
448 sprintf (pout," Minute = %i -- greater than 59\n",val);
449 gaprnt (0,pout);
450 return (NULL);
451 }
452 if (tolower(*ch)!='z') {
453 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
454 gaprnt (0," 'z' delimiter is missing \n");
455 return (NULL);
456 }
457 dtim->mn = val;
458 ch++;
459 if (*ch>='0' && *ch<='9') ch = intprs (ch,&val);
460 else val = def->dy;
461 } else {
462 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
463 gaprnt (0," Missing minute value \n");
464 return (NULL);
465 }
466 } else {
467 ch++;
468 if (*ch>='0' && *ch<='9') ch = intprs (ch,&val);
469 else val = def->dy;
470 }
471 } else flag = 2;
472 dtim->dy = val;
473 } else flag = 1;
474
475 monam[0] = tolower(*ch);
476 monam[1] = tolower(*(ch+1));
477 monam[2] = tolower(*(ch+2));
478 monam[3] = '\0';
479
480 i = 0;
481 while (i<12 && !cmpwrd(monam,mons[i]) ) i++;
482 i++;
483
484 if (i==13) {
485 if (flag==1) {
486 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
487 gaprnt (0," Expected month abbreviation, none found\n");
488 return (NULL);
489 }
490 if (flag==2) {
491 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
492 gaprnt (0," Missing month abbreviation or 'z' delimiter\n");
493 return (NULL);
494 }
495 dtim->mo = def->mo;
496 dtim->yr = def->yr;
497 } else {
498 dtim->mo = i;
499 ch+=3;
500 /*mf ---- parse year ---- mf */
501 if (*ch>='0' && *ch<='9') {
502 /*mf 971117 --- use fullyear only if year 1 = 0001 mf*/
503 if(*(ch+2)>='0' && *(ch+2)<='9') {
504 mfcmn.fullyear=1;
505 } else {
506 mfcmn.fullyear=0;
507 }
508 ch = intprs (ch,&val);
509 } else {
510 val = def->yr;
511 }
512
513 /*mf --- turn off setting of < 100 years to 1900 or 2000 ---mf*/
514 if(mfcmn.fullyear == 0) {
515 if (val<50) val+=2000;
516 else if (val<100) val+=1900;
517 }
518 dtim->yr = val;
519 }
520
521 i = mosiz[dtim->mo];
522 if (dtim->mo==2 && qleap(dtim->yr)) i = 29;
523 if (dtim->dy > i) {
524 gaprnt (0,"Syntax Error: Invalid Date/Time value.\n");
525 sprintf (pout," Day = %i -- greater than %i \n",dtim->dy,i);
526 gaprnt (0,pout);
527 return (NULL);
528 }
529 return (ch);
530 }
531
532 /* Parse a relative date/time (offset). Format is:
533
534 nn (yr/mo/dy/hr/mn)
535
536 Examples: 5mo
537 1dy12hr
538 etc.
539
540 Missing values are filled in with 0s. */
541
rdtprs(char * ch,struct dt * dtim)542 char *rdtprs (char *ch, struct dt *dtim) {
543 int flag,val;
544 char *pos;
545 char id[3];
546
547 pos = ch;
548
549 dtim->yr = 0;
550 dtim->mo = 0;
551 dtim->dy = 0;
552 dtim->hr = 0;
553 dtim->mn = 0;
554
555 flag = 1;
556
557 while (*ch>='0' && *ch<='9') {
558 flag = 0;
559 ch = intprs(ch,&val);
560 id[0] = *ch; id[1] = *(ch+1); id[2] = '\0';
561 if (cmpwrd("yr",id)) dtim->yr = val;
562 else if (cmpwrd("mo",id)) dtim->mo = val;
563 else if (cmpwrd("dy",id)) dtim->dy = val;
564 else if (cmpwrd("hr",id)) dtim->hr = val;
565 else if (cmpwrd("mn",id)) dtim->mn = val;
566 else {
567 gaprnt (0,"Syntax Error: Invalid Date/Time offset.\n");
568 sprintf (pout," Expecting yr/mo/dy/hr/mn, found %s\n",id);
569 gaprnt (0,pout);
570 return (NULL);
571 }
572 ch+=2;
573 }
574 if (flag) {
575 gaprnt (0,"Syntax Error: Invalid Date/Time offset.\n");
576 gaprnt (0," No offset value given\n");
577 return (NULL);
578 }
579 return (ch);
580 }
581 static char num[10]={'0','1','2','3','4','5','6','7','8','9'};
582
gaedit(float val,char * chars,int icode)583 int gaedit (float val, char *chars, int icode) {
584 int i,j,k,sign;
585 char ch[20], *cc;
586 float vvv;
587 int w1,w2;
588
589 vvv=icode;
590 vvv=pow(10.0,vvv);
591 vvv=vvv*val;
592
593 sign=0;
594 if (vvv<0.0) {sign=1; vvv=vvv*(-1.0);}
595
596 i = 1;
597 ch[19] = '\0';
598 cc = &ch[18];
599 w2 = vvv+0.5;
600 while (w2>0) {
601 w1 = w2/10;
602 w1 = w1*10;
603 w1 = w2-w1;
604 *cc = num[w1];
605 cc--; i++;
606 icode--;
607 k=0;
608 if (icode==0) {
609 *cc='.';
610 cc--;
611 i++;
612 k=1;
613 }
614 w2 = w2/10;
615 }
616 if (k) {i++; *cc='0'; cc--;}
617 if (i==1) {i++; *cc='0';}
618 else {
619 if (sign) {*cc='-'; i++;}
620 else *cc++;
621 }
622 for (j=0;j<i;j++) {
623 *chars = *cc;
624 cc++; chars++;
625 }
626 return (i-1);
627 }
628
629 /* Compares two strings. A match occurs if the leading
630 blank-delimited words in the two strings match. CR and NULL also
631 serve as delimiters. */
632
cmpwrd(char * ch1,char * ch2)633 int cmpwrd (char *ch1, char *ch2) {
634
635 while (*ch1==' '||*ch1=='\t') ch1++; /* Advance past leading blanks. */
636 while (*ch2==' '||*ch2=='\t') ch2++;
637
638 while (*ch1 == *ch2) {
639 if (*ch1==' '||*ch1=='\t'||*ch1=='\0'||*ch1=='\n'||*ch1=='\r' ) return (1);
640 ch1++; ch2++;
641 }
642
643 if ( (*ch1==' '||*ch1=='\t'||*ch1=='\0'||*ch1=='\n'||*ch1=='\r') &&
644 (*ch2==' '||*ch2=='\t'||*ch2=='\0'||*ch2=='\n'||*ch2=='\r') ) return (1);
645 return (0);
646 }
647 /* case insensitive version of cmpwrd */
648
cmpwrdl(char * ch1,char * ch2)649 int cmpwrdl (char *ch1, char *ch2) {
650 if(ch1 == NULL || ch2 == NULL) return(0);
651
652 while (*ch1==' '||*ch1=='\t') ch1++; /* Advance past leading blanks. */
653 while (*ch2==' '||*ch2=='\t') ch2++;
654
655 while (tolower(*ch1) == tolower(*ch2)) {
656 if (*ch1==' '||*ch1=='\t'||*ch1=='\0'||*ch1=='\n'||*ch1=='\r' ) return (1);
657 ch1++; ch2++;
658 }
659
660 if ( (*ch1==' '||*ch1=='\t'||*ch1=='\0'||*ch1=='\n'||*ch1=='\r' ) &&
661 (*ch2==' '||*ch2=='\t'||*ch2=='\0'||*ch2=='\n'||*ch2=='\r' ) ) return (1);
662 return (0);
663 }
664
665 /* Moves a pointer to the start of the next blank-delimited word
666 in a string. If not found, NULL is returned. */
667
nxtwrd(char * ch)668 char * nxtwrd (char *ch) {
669
670 while (*ch!=' '&&*ch!='\t') { /* Skip 1st word */
671 if (*ch == '\0' || *ch == '\n' || *ch == '\r') return (NULL);
672 ch++;
673 }
674 while (*ch==' '||*ch=='\t') ch++; /* Find next word */
675 if (*ch == '\0' || *ch == '\n' || *ch == '\r') return (NULL);
676 return (ch);
677 }
678
679 /* Figures out how many digits are after the decimal point */
680
gxdeci(float val)681 int gxdeci (float val) {
682 int code,ival;
683
684 code = 0;
685 val = fabs(val);
686 ival = val;
687 val = val - (float)ival;
688 while (code<4&&val>0.01) {
689 val=val*10.0;
690 ival=val;
691 val=val-(float)ival;
692 code++;
693 }
694 return (code);
695 }
696
697 /* Linear conversion routine for dimension conversions. */
698
liconv(float * vals,float v)699 float liconv (float *vals, float v) {
700 return ( (*vals * v) + *(vals+1) );
701 }
702
703 /* Non-linear scaling routine for discrete levels. Linear interp
704 between levels. Scaling beyond upper and lower bounds is
705 linear based on the last and first grid spacing, respectively.
706 In each case a pointer to a list of values is provided. The
707 list contains in its first element the number of values
708 in the list. */
709
710 /* Convert a grid value to a world coordinate value.
711 This operation needs to be efficient, since it gets done
712 very often. */
713
gr2lev(float * vals,float gr)714 float gr2lev (float *vals, float gr) {
715 int i, num;
716 float x;
717 if (gr<1.0) return ( *(vals+1) + (1.0-gr)*(*(vals+1)-*(vals+2)) );
718 if (gr>*vals) {
719 i = (int)(*vals+0.1);
720 return ( *(vals+i) + (gr-*vals)*(*(vals+i)-*(vals+i-1)) );
721 }
722 i = (int)gr;
723 return (*(vals+i)+((gr-(float)i)*(*(vals+i+1)-*(vals+i))));
724 }
725
726 /* Convert from world coordinate value to grid value. This operation
727 is not set up to be efficient, under the assumption that it won't
728 get done all that often. */
729
lev2gr(float * vals,float lev)730 float lev2gr (float *vals, float lev) {
731 int i,num;
732 float gr;
733 num = (int)(*vals+0.1);
734 for (i=1; i<num; i++) {
735 if ( (lev >= *(vals+i) && lev <= *(vals+i+1)) ||
736 (lev <= *(vals+i) && lev >= *(vals+i+1)) ) {
737 gr = (float)i + (lev - *(vals+i))/(*(vals+i+1) - *(vals+i));
738 return (gr);
739 }
740 }
741 if (*(vals+1)<*(vals+num)) {
742 if (lev<*(vals+1)) {
743 gr = 1.0 + ((lev-*(vals+1))/(*(vals+2)-*(vals+1)));
744 return (gr);
745 }
746 gr = (float)i + ((lev-*(vals+i))/(*(vals+i)-*(vals+i-1)));
747 return (gr);
748 } else {
749 if (lev>*(vals+1)) {
750 gr = 1.0 + ((lev-*(vals+1))/(*(vals+2)-*(vals+1)));
751 return (gr);
752 }
753 gr = (float)i + ((lev-*(vals+i))/(*(vals+i)-*(vals+i-1)));
754 return (gr);
755 }
756 }
757
758 /* Parses a number in a character string.
759 This routine will detect numbers of the form:
760 nnnn
761 -nnnn
762
763 Args: ch - pointer to the number, in character form.
764 val - integer value returned
765 return value - address of 1st character past the
766 number parsed. NULL if no number found
767 at pointer ch or if the number is an
768 invalid format.
769 */
770
intprs(char * ch,int * val)771 char *intprs (char *ch, int *val) {
772
773 int nflag,flag;
774
775 nflag = 0;
776 if (*ch=='-') { nflag = 1; ch++; }
777 else if (*ch=='+') ch++;
778
779 *val = 0;
780 flag = 1;
781
782 while (*ch>='0' && *ch<='9') {
783 *val = *val*10 + (int)(*ch-'0');
784 flag = 0;
785 ch++;
786 }
787
788 if (flag) return (NULL);
789
790 if (nflag) *val = -1 * *val;
791 return (ch);
792 }
793
longprs(char * ch,long * val)794 char *longprs (char *ch, long *val) {
795
796 int nflag,flag;
797
798 nflag = 0;
799 if (*ch=='-') { nflag = 1; ch++; }
800 else if (*ch=='+') ch++;
801
802 *val = 0;
803 flag = 1;
804
805 while (*ch>='0' && *ch<='9') {
806 *val = *val*10 + (int)(*ch-'0');
807 flag = 0;
808 ch++;
809 }
810
811 if (flag) return (NULL);
812
813 if (nflag) *val = -1 * *val;
814 return (ch);
815 }
816
817
818 /* Parses a number in a character string.
819 This routine will detect numbers of the form:
820 nnnn
821 nnnn.fff
822 nnnn.fffExxx
823
824 Args: ch - pointer to the number, in character form.
825 val - floating point value returned
826 return value - address of 1st character past the
827 number parsed. NULL if no number found
828 at pointer ch or if the number is an
829 invalid format.
830 */
831
832 #ifdef HAVE_STRTOD
833
valprs(char * ch,float * val)834 char *valprs (char *ch, float *val) {
835 char * pos;
836
837 *val = (float)strtod(ch, &pos);
838 if (pos==ch) {
839 return NULL;
840 } else {
841 return pos;
842 }
843 }
844
845 #else
846
valprs(char * ch,float * val)847 char *valprs (char *ch, float *val) {
848
849 int nflag,dflag,eflag,enflag,flag,cont;
850 int pflag,epflag,evflag;
851 float exp,dfp;
852 int zip;
853
854 flag=0;
855 nflag=0;dflag=0;eflag=0;enflag=0;
856 pflag=0;epflag=0;evflag=0;
857 *val=0.0;exp=0.0;dfp=0.1;
858 zip='0';
859
860 cont=1;
861 while (cont) {
862
863 if (*ch>='0' && *ch<='9') {
864 if (!flag) flag=1;
865 if (eflag) {evflag=1; exp=(exp*10)+(*ch-zip);}
866 else if (dflag) { *val = *val+((*ch-zip)*dfp); dfp=dfp/10.0; }
867 else *val = (*val*10.0)+(*ch-zip);
868
869 } else if (*ch=='-') {
870 if (eflag&&!evflag) {
871 if (enflag) {cont=0; flag=0;}
872 enflag=1;
873 }
874 else if (!flag) {
875 if (nflag) {cont=0; flag=0;}
876 nflag=1;
877 } else cont=0;
878
879 } else if (*ch=='+') {
880 if (eflag&&!evflag) {
881 if (epflag) {cont=0; flag=0;}
882 epflag=1;
883 }
884 else if (!flag) {
885 if (pflag) {cont=0; flag=0;}
886 pflag=1;
887 } else cont=0;
888
889 } else if (*ch=='.') {
890 if (dflag||eflag) {cont=0;}
891 else dflag=1;
892
893 } else if (*ch=='e' || *ch=='E') {
894 if (eflag) {
895 cont=0;
896 if (!evflag) flag=0;
897 }
898 else if (flag) {eflag=1; dflag=0;}
899 else cont=0;
900
901 } else cont=0;
902
903 if (cont) ch++;
904 }
905
906 if (flag) {
907 if (nflag) *val = *val*(-1.0);
908 if (eflag) {
909 if (enflag) exp = exp*(-1.0);
910 *val = *val*(pow(10.0,exp));
911 }
912 return (ch);
913 } else return (NULL);
914
915 }
916
917 #endif
918
919
920 /* dimprs parses a dimension 'expression', ie, where the user
921 specifies an absolute or relative dimension value.
922 The format is:
923
924 dim op val
925
926 where: dim = x,y,z,t,lat,lon,lev,time
927 op = +, -, or =
928 val = dimension value
929
930 Examples are:
931 t=1
932 lev=500
933 time=00:00z01jan1989
934
935 The coordinates are evaluated with respect to the coordinate
936 transformations of the file descriptor passed. Grid
937 coordinates are returned. Relative offsets are evaluated
938 from the values in the status block.
939
940 In addition, r=radius is also supported. The dimension value
941 returned is the radius, and the dimension number returned
942 is 9. This is only valid for stn type files.
943
944 wflag is set to zero if the dimension expression was
945 grid coordinates; 1 if it was world coordinates.
946 */
947
dimprs(char * ch,struct gastat * pst,struct gafile * pfi,int * dim,float * d,int type,int * wflag)948 char *dimprs (char *ch, struct gastat *pst, struct gafile *pfi,
949 int *dim, float *d, int type, int *wflag) {
950 char *pos, cc, cc2, cd, *frst;
951 char name[5];
952 int dval, flag,i,op;
953 struct dt dtim;
954 float (*conv) (float *, float);
955 float *cvals,v;
956
957 frst = ch;
958 i = 0;
959 while (*ch>='a' && *ch<='z' && i<6) {
960 name[i] = *ch;
961 ch++; i++;
962 }
963 name[i] = '\0';
964 if (i>4) {
965 gaprnt (0,"Syntax Error: Invalid dimension expression\n");
966 sprintf (pout," Expecting X/Y/Z/T/LAT/LON/LEV/TIME, found %s\n",
967 name);
968 gaprnt (0,pout);
969 return (NULL);
970 }
971
972 if (*ch == '=') op = 0;
973 else if (*ch == '+') op = 1;
974 else if (*ch == '-') op = 2;
975 else {
976 gaprnt (0,"Syntax Error: Invalid dimension expression\n");
977 sprintf (pout," Expecting +/-/= operator, found %c\n",*ch);
978 gaprnt (0,pout);
979 return (NULL);
980 }
981
982 ch++;
983 if (cmpwrd("time",name)) {
984 if (op==0) {
985 if ( (pos=adtprs(ch,&(pst->tmin),&dtim)) == NULL ) {
986 gaprnt (0," Invalid dimension expression\n");
987 return (NULL);
988 }
989 } else {
990 if ( (pos=rdtprs(ch,&dtim)) == NULL ) {
991 return (NULL);
992 }
993 }
994 } else {
995 if ( (pos=valprs(ch,&v)) == NULL ) {
996 gaprnt (0,"Syntax Error: Invalid dimension expression\n");
997 gaprnt (0," Dimension value missing or invalid\n");
998 return (NULL);
999 }
1000 }
1001 ch = pos;
1002
1003 /* We now have all the info we need about this dimension
1004 expression. We can now evaluate it. */
1005
1006 if (cmpwrd("x",name)) *dim = 0;
1007 else if (cmpwrd("y",name)) *dim = 1;
1008 else if (cmpwrd("z",name)) *dim = 2;
1009 else if (cmpwrd("t",name)) *dim = 3;
1010 else if (cmpwrd("lon",name)) *dim = 4;
1011 else if (cmpwrd("lat",name)) *dim = 5;
1012 else if (cmpwrd("lev",name)) *dim = 6;
1013 else if (cmpwrd("time",name)) *dim = 7;
1014 else if (type==0 && cmpwrd("r",name)) *dim = 9;
1015 else {
1016 gaprnt (0,"Syntax Error: Invalid dimension expression\n");
1017 sprintf (pout," Expecting X/Y/Z/T/LAT/LON/LEV/TIME, found %s\n",
1018 name);
1019 gaprnt (0,pout);
1020 return (NULL);
1021 }
1022
1023 *wflag = 0;
1024 if (*dim == 9) {
1025 *d = v;
1026 return (ch);
1027 }
1028
1029 if (*dim < 4) {
1030 if (op==0) {
1031 *d = v + pfi->dimoff[*dim];
1032 return (ch);
1033 } else {
1034 if (*dim==pst->idim || *dim==pst->jdim) {
1035 gaprnt (0,"Syntax Error: Invalid dimension expression\n");
1036 gaprnt (0," Cannot use an offset value with a varying ");
1037 gaprnt (0,"dimension\n");
1038 sprintf (pout," Varying dimension = %i \n",*dim);
1039 gaprnt (0,pout);
1040 return (NULL);
1041 }
1042 if (*dim < 3) {
1043 if (pfi->type==1) {
1044 conv = pfi->ab2gr[*dim];
1045 cvals = pfi->abvals[*dim];
1046 *d = conv(cvals,pst->dmin[*dim]);
1047 } else {
1048 *d = pst->dmin[*dim];
1049 }
1050 } else {
1051 *d = t2gr(pfi->abvals[3],&(pst->tmin));
1052 }
1053 if (op==1) *d = *d + v;
1054 if (op==2) *d = *d - v;
1055 return (ch);
1056 }
1057 } else {
1058 *dim = *dim - 4;
1059 *wflag = 1;
1060 if (op>0) {
1061 if (*dim==pst->idim || *dim==pst->jdim) {
1062 gaprnt (0,"Syntax Error: Invalid dimension expression\n");
1063 gaprnt (0," Cannot use an offset value with a varying ");
1064 gaprnt (0,"dimension\n");
1065 sprintf (pout," Varying dimension = %i \n",*dim);
1066 gaprnt (0,pout);
1067 return (NULL);
1068 }
1069 if (*dim==3) {
1070 if (op==1) timadd (&(pst->tmin),&dtim);
1071 if (op==2) timsub (&(pst->tmin),&dtim);
1072 } else {
1073 if (op==1) v = pst->dmin[*dim] + v;
1074 if (op==2) v = pst->dmin[*dim] - v;
1075 }
1076 }
1077 if (*dim < 3) {
1078 if (pfi->type==1 || pfi->type==4) {
1079 conv = pfi->ab2gr[*dim];
1080 cvals = pfi->abvals[*dim];
1081 *d = conv(cvals,v);
1082 } else {
1083 *d = v;
1084 }
1085 } else {
1086 *d = t2gr(pfi->abvals[3],&dtim);
1087 }
1088 return (ch);
1089 }
1090 }
1091
1092 /*mf version
1093 convert all upper case alphabetic characters to lower case.
1094 The GrADS system is case insensitive, and assumes lower case
1095 internally in most cases. */
1096
lowcas(char * ch)1097 void lowcas (char *ch) {
1098 int i;
1099 int qflag=0;
1100
1101 while (*ch!='\0' && *ch!='\n') {
1102 i = *ch;
1103 /*
1104 mf - 940415 : do not turn to lower case if in "'s
1105 */
1106 if(*ch == '\"' && qflag == 0 ) {
1107 qflag=1;
1108 } else if(*ch == '\"' && qflag == 1 ) {
1109 qflag=0;
1110 }
1111 if (i>64 && i<91 && qflag==0) {
1112 i+=32;
1113 *ch = i;
1114 } else if(i == 95) {
1115 *ch=i;
1116 }
1117 ch++;
1118 }
1119 }
1120
1121
1122 /* convert all upper case alphabetic characters to lower case.
1123 The GrADS system is case insensitive, and assumes lower case
1124 internally in most cases. */
1125 /*----------original version
1126 void lowcas (char *ch) {
1127 int i;
1128
1129 while (*ch!='\0' && *ch!='\n') {
1130 i = *ch;
1131 if (i>64 && i<91) {
1132 i+=32;
1133 *ch = i;
1134 }
1135 ch++;
1136 }
1137 }
1138 ------------original version */
1139
1140 /* convert to upper case */
1141
uppcas(char * ch)1142 void uppcas (char *ch) {
1143 int i;
1144
1145 while (*ch!='\0' && *ch!='\n') {
1146 i = *ch;
1147 if (i>96 && i<123) {
1148 i -= 32;
1149 *ch = i;
1150 }
1151 ch++;
1152 }
1153 }
1154
1155 /* Copies a string of a specified length, or when \0 or \n is hit.
1156 Trailing blanks are removed, and the output string is terminated
1157 with '\0'. */
1158
getstr(char * ch1,char * ch2,int len)1159 void getstr (char *ch1, char *ch2, int len) {
1160 char *ch;
1161
1162 ch = ch1;
1163 while (len>0 && *ch2!='\n' && *ch2!='\0') {
1164 *ch1 = *ch2;
1165 len--;
1166 ch1++; ch2++;
1167 }
1168 ch1--;
1169 while (ch1>=ch && *ch1==' ') ch1--;
1170 ch1++;
1171 *ch1 = '\0';
1172 }
1173
1174 /* Copies a word of a specified length, or when \0 or \n or \r or ' ' is
1175 encountered. The word is terminated with '\0'. */
1176
getwrd(char * ch1,char * ch2,int len)1177 void getwrd (char *ch1, char *ch2, int len) {
1178 char *ch;
1179
1180 ch = ch1;
1181 while (len>0 && *ch2!='\n' && *ch2!='\0' && *ch2!='\r' && *ch2!=' ' ) {
1182 *ch1 = *ch2;
1183 len--;
1184 ch1++; ch2++;
1185 }
1186 *ch1 = '\0';
1187 }
1188
1189 /* Determines word length up to next delimiter */
1190
wrdlen(char * ch2)1191 int wrdlen (char *ch2) {
1192 int len;
1193 len = 0;
1194 while (*ch2!='\n' && *ch2!='\0' && *ch2!=' ' && *ch2!='\t') {
1195 len++;
1196 ch2++;
1197 }
1198 return(len);
1199 }
1200
1201 /* Get minimum and maximum grid value. Set rmin and rmax in the
1202 grid descriptor. */
1203
gamnmx(struct gagrid * pgr)1204 void gamnmx (struct gagrid *pgr) {
1205 int i,size;
1206 float *r;
1207 int cnt;
1208 size = pgr->isiz * pgr->jsiz;
1209 if (size==1) return;
1210 pgr->rmin= 9.99E35;
1211 pgr->rmax= -9.99E35;
1212 r = pgr->grid;
1213 cnt=0;
1214 for (i=0;i<size;i++) {
1215 if (*r != pgr->undef) {
1216 cnt++;
1217 if (pgr->rmin>*r) pgr->rmin = *r;
1218 if (pgr->rmax<*r) pgr->rmax = *r;
1219 }
1220 r++;
1221 }
1222 if (cnt==0 || pgr->rmin==9.99e35 || pgr->rmax==-9.99e35) {
1223 pgr->rmin = pgr->undef;
1224 pgr->rmax = pgr->undef;
1225 }
1226 }
1227
1228 /* Determine min and max values of station data */
1229
gasmnmx(struct gastn * stn)1230 void gasmnmx (struct gastn *stn) {
1231 struct garpt *rpt;
1232
1233 stn->smin = stn->undef;
1234 stn->smax = stn->undef;
1235 rpt = stn->rpt;
1236 while (rpt!=NULL) {
1237 if (rpt->val != stn->undef) {
1238 if (stn->smin == stn->undef) stn->smin = rpt->val;
1239 else {
1240 if (stn->smin > rpt->val) stn->smin = rpt->val;
1241 }
1242 if (stn->smax == stn->undef) stn->smax = rpt->val;
1243 else {
1244 if (stn->smax < rpt->val) stn->smax = rpt->val;
1245 }
1246 }
1247 rpt = rpt->rpt;
1248 }
1249 }
1250
1251 /*mf version Remove blanks from a string */
garemb(char * ch)1252 int garemb (char *ch) {
1253 char *cc;
1254 int cnt;
1255 int qflag=0;
1256
1257 cc = ch;
1258 cnt = 0;
1259
1260 while ( *ch!='\n' && *ch!='\0' ) {
1261
1262 /*
1263 mf - 940415 : do not remove blanks if string in "'s
1264 */
1265 if(*ch == '\"' && qflag == 0 ) {
1266 qflag=1;
1267 } else if(*ch == '\"' && qflag == 1 ) {
1268 qflag=0;
1269 }
1270
1271 if ( ( (*ch!=' ') || qflag ) && (*ch!='\"') ) {
1272 *cc = *ch;
1273 cc++; cnt++;
1274 }
1275 ch++;
1276 }
1277 *cc = '\0';
1278 return (cnt);
1279 }
1280
1281 /*---- original ------
1282
1283 int garemb (char *ch) {
1284 char *cc;
1285 int cnt;
1286
1287 cc = ch;
1288 cnt = 0;
1289 while ( *ch!='\n' && *ch!='\0' ) {
1290 if (*ch!=' ') {
1291 *cc = *ch;
1292 cc++; cnt++;
1293 }
1294 ch++;
1295 }
1296 *cc = '\0';
1297 return (cnt);
1298 }
1299 */
1300
1301 static float glts15[40] = {
1302 -86.60,-82.19,-77.76,-73.32,-68.88,-64.43,-59.99,
1303 -55.55,-51.11,-46.66,-42.22,-37.77,-33.33,-28.89,
1304 -24.44,-20.00,-15.55,-11.11, -6.67, -2.22, 2.22,
1305 6.67, 11.11, 15.55, 20.00, 24.44, 28.89, 33.33,
1306 37.77, 42.22, 46.66, 51.11, 55.55, 59.99, 64.43,
1307 68.88, 73.32, 77.76, 82.19, 86.60};
1308 static float glts20[52] = {
1309 -87.38,-83.98,-80.56,-77.13,-73.71,-70.28,-66.85,
1310 -63.42,-59.99,-56.57,-53.14,-49.71,-46.28,-42.85,
1311 -39.43,-36.00,-32.57,-29.14,-25.71,-22.28,-18.86,
1312 -15.43,-12.00, -8.57, -5.14, -1.71, 1.71, 5.14,
1313 8.57, 12.00, 15.43, 18.86, 22.28, 25.71, 29.14,
1314 32.57, 36.00, 39.43, 42.85, 46.28, 49.71, 53.14,
1315 56.57, 59.99, 63.42, 66.85, 70.28, 73.71, 77.13,
1316 80.56, 83.98, 87.38};
1317 static float glts30[80] = {
1318 -88.29, -86.07, -83.84, -81.61, -79.37, -77.14, -74.90,
1319 -72.67, -70.43, -68.20, -65.96, -63.72, -61.49, -59.25,
1320 -57.02, -54.78, -52.55, -50.31, -48.07, -45.84, -43.60,
1321 -41.37, -39.13, -36.89, -34.66, -32.42, -30.19, -27.95,
1322 -25.71, -23.48, -21.24, -19.01, -16.77, -14.53, -12.30,
1323 -10.06, -7.83, -5.59, -3.35, -1.12, 1.12, 3.35,
1324 5.59, 7.83, 10.06, 12.30, 14.53, 16.77, 19.01,
1325 21.24, 23.48, 25.71, 27.95, 30.19, 32.42, 34.66,
1326 36.89, 39.13, 41.37, 43.60, 45.84, 48.07, 50.31,
1327 52.55, 54.78, 57.02, 59.25, 61.49, 63.72, 65.96,
1328 68.20, 70.43, 72.67, 74.90, 77.14, 79.37, 81.61,
1329 83.84, 86.07, 88.29};
1330
1331 static float glats[102] = {
1332 -88.66,-86.91,-85.16,-83.41,-81.65,-79.90,-78.14,-76.39,-74.63,
1333 -72.88,-71.12,-69.36,-67.61,-65.85,-64.10,-62.34,-60.58,-58.83,
1334 -57.07,-55.32,-53.56,-51.80,-50.05,-48.29,-46.54,-44.78,-43.02,
1335 -41.27,-39.51,-37.76,-36.00,-34.24,-32.49,-30.73,-28.98,-27.22,
1336 -25.46,-23.71,-21.95,-20.19,-18.44,-16.68,-14.93,-13.17,-11.41,
1337 -9.66, -7.90, -6.15, -4.39, -2.63, -0.88, 0.88, 2.63, 4.39,
1338 6.15, 7.90, 9.66, 11.41, 13.17, 14.93, 16.68, 18.44, 20.19,
1339 21.95, 23.71, 25.46, 27.22, 28.98, 30.73, 32.49, 34.24, 36.00,
1340 37.76, 39.51, 41.27, 43.02, 44.78, 46.54, 48.29, 50.05, 51.80,
1341 53.56, 55.32, 57.07, 58.83, 60.58, 62.34, 64.10, 65.85, 67.61,
1342 69.36, 71.12, 72.88, 74.63, 76.39, 78.14, 79.90, 81.65, 83.41,
1343 85.16, 86.91, 88.66 };
1344
1345 static float m32lts[32] = {-20.453, -18.01, -15.763, -13.738,
1346 -11.95, -10.405, -9.097, -8.010, -7.120, -6.392, -5.253, -4.25,
1347 -3.25, -2.25, -1.25, -0.25, 0.25, 1.25, 2.25, 3.25, 4.25, 5.253,
1348 6.392, 7.12, 8.01, 9.097, 10.405, 11.95, 13.738, 15.763, 18.01,
1349 20.453};
1350
1351 /* From Mike Timlin */
1352 static float gltst62[94] = {
1353 -88.542, -86.6531, -84.7532, -82.8508, -80.9473, -79.0435,
1354 -77.1394, -75.2351, -73.3307, -71.4262, -69.5217, -67.6171,
1355 -65.7125, -63.8079, -61.9033, -59.9986, -58.0939, -56.1893,
1356 -54.2846, -52.3799, -50.4752, -48.5705, -46.6658, -44.7611,
1357 -42.8564, -40.9517, -39.047, -37.1422, -35.2375, -33.3328,
1358 -31.4281, -29.5234, -27.6186, -25.7139, -23.8092, -21.9044,
1359 -19.9997, -18.095, -16.1902, -14.2855, -12.3808, -10.47604,
1360 -8.57131, -6.66657, -4.76184, -2.8571, -0.952368, 0.952368,
1361 2.8571, 4.76184, 6.66657, 8.57131, 10.47604, 12.3808, 14.2855,
1362 16.1902, 18.095, 19.9997, 21.9044, 23.8092, 25.7139, 27.6186,
1363 29.5234, 31.4281, 33.3328, 35.2375, 37.1422, 39.047, 40.9517,
1364 42.8564, 44.7611, 46.6658, 48.5705, 50.4752, 52.3799, 54.2846,
1365 56.1893, 58.0939, 59.9986, 61.9033, 63.8079, 65.7125, 67.6171,
1366 69.5217, 71.4262, 73.3307, 75.2351, 77.1394, 79.0435, 80.9473,
1367 82.8508, 84.7532, 86.6531, 88.542 };
1368
1369 /* Given the starting point and the length, return the MOM32 lats */
1370
gamo32(int istrt,int num)1371 float *gamo32 (int istrt, int num) {
1372 int size;
1373 float *vals;
1374
1375 istrt--;
1376 if (istrt+num > 32) {
1377 gaprnt (0,"Open Error: Invalid MOM32 scaling.\n");
1378 gaprnt (0," Maximum 32 latitudes exceeded \n");
1379 return (NULL);
1380 }
1381 size = (num+3) * sizeof(float);
1382 vals = (float *)malloc(size);
1383 if (vals==NULL) {
1384 gaprnt (0,"Memory Allocation Error: MOM32 Grid Scaling\n");
1385 return (NULL);
1386 }
1387 *vals = (float)num;
1388 for (size=0; size<num; size++) *(vals+size+1) = m32lts[size+istrt];
1389 *(vals+num+1) = -999.9;
1390 return (vals);
1391 }
1392
1393
1394 /* Given the starting point and the length, return the gaussian lats */
1395
gagaus(int istrt,int num)1396 float *gagaus (int istrt, int num) {
1397 int size;
1398 float *vals;
1399
1400 istrt--;
1401 if (istrt+num > 102) {
1402 gaprnt (0,"Open Error: Invalid GAUSR40 scaling.\n");
1403 gaprnt (0," Maximum 102 latitudes exceeded \n");
1404 return (NULL);
1405 }
1406 size = (num+3) * sizeof(float);
1407 vals = (float *)malloc(size);
1408 if (vals==NULL) {
1409 gaprnt (0,"Memory Allocation Error: Gaussian Grid Scaling\n");
1410 return (NULL);
1411 }
1412 *vals = (float)num;
1413 for (size=0; size<num; size++) *(vals+size+1) = glats[size+istrt];
1414 *(vals+num+1) = -999.9;
1415 return (vals);
1416 }
1417
1418 /* Given the starting point and the length, return the gaussian lats
1419 for R30 grids */
1420
gags30(int istrt,int num)1421 float *gags30 (int istrt, int num) {
1422 int size;
1423 float *vals;
1424
1425 istrt--;
1426 if (istrt+num > 80) {
1427 gaprnt (0,"Open Error: Invalid GAUSR30 scaling.\n");
1428 gaprnt (0," Maximum 80 latitudes exceeded \n");
1429 return (NULL);
1430 }
1431 size = (num+3) * sizeof(float);
1432 vals = (float *)malloc(size);
1433 if (vals==NULL) {
1434 gaprnt (0,"Memory Allocation Error: Gaussian Grid Scaling\n");
1435 return (NULL);
1436 }
1437 *vals = (float)num;
1438 for (size=0; size<num; size++) *(vals+size+1) = glts30[size+istrt];
1439 *(vals+num+1) = -999.9;
1440 return (vals);
1441 }
1442
1443 /* Given the starting point and the length, return the gaussian lats
1444 for R20 grids */
1445
gags20(int istrt,int num)1446 float *gags20 (int istrt, int num) {
1447 int size;
1448 float *vals;
1449
1450 istrt--;
1451 if (istrt+num > 52) {
1452 gaprnt (0,"Open Error: Invalid GAUSR20 scaling.\n");
1453 gaprnt (0," Maximum 52 latitudes exceeded \n");
1454 return (NULL);
1455 }
1456 size = (num+3) * sizeof(float);
1457 vals = (float *)malloc(size);
1458 if (vals==NULL) {
1459 gaprnt (0,"Memory Allocation Error: Gaussian Grid Scaling\n");
1460 return (NULL);
1461 }
1462 *vals = (float)num;
1463 for (size=0; size<num; size++) *(vals+size+1) = glts20[size+istrt];
1464 *(vals+num+1) = -999.9;
1465 return (vals);
1466 }
1467
1468 /* Given the starting point and the length, return the gaussian lats
1469 for R15 grids */
1470
gags15(int istrt,int num)1471 float *gags15 (int istrt, int num) {
1472 int size;
1473 float *vals;
1474
1475 istrt--;
1476 if (istrt+num > 40) {
1477 gaprnt (0,"Open Error: Invalid GAUSR15 scaling.\n");
1478 gaprnt (0," Maximum 40 latitudes exceeded \n");
1479 return (NULL);
1480 }
1481 size = (num+3) * sizeof(float);
1482 vals = (float *)malloc(size);
1483 if (vals==NULL) {
1484 gaprnt (0,"Memory Allocation Error: Gaussian Grid Scaling\n");
1485 return (NULL);
1486 }
1487 *vals = (float)num;
1488 for (size=0; size<num; size++) *(vals+size+1) = glts15[size+istrt];
1489 *(vals+num+1) = -999.9;
1490 return (vals);
1491 }
1492
1493 /* Given the starting point and the length, return the gaussian lats
1494 for T62 grids */
1495 /* From Mike Timlin */
1496
gagst62(int istrt,int num)1497 float *gagst62 (int istrt, int num) {
1498 int size;
1499 float *vals;
1500
1501 istrt--;
1502 if (istrt+num > 94) {
1503 gaprnt (0,"Open Error: Invalid GAUST62 scaling.\n");
1504 gaprnt (0," Maximum 94 latitudes exceeded \n");
1505 return (NULL);
1506 }
1507 size = (num+3) * sizeof(float);
1508 vals = (float *)malloc(size);
1509 if (vals==NULL) {
1510 gaprnt (0,"Memory Allocation Error: Gaussian Grid Scaling\n");
1511 return (NULL);
1512 }
1513 *vals = (float)num;
1514 for (size=0; size<num; size++) *(vals+size+1) = gltst62[size+istrt];
1515 *(vals+num+1) = -999.9;
1516 return (vals);
1517 }
1518
1519 char *monc[12] = {"JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG",
1520 "SEP","OCT","NOV","DEC"};
1521
gat2ch(struct dt * dtim,int tinc,char * ch)1522 int gat2ch (struct dt *dtim, int tinc, char *ch) {
1523 int mn1,mn2,hr1,hr2,dy1,dy2,len,mnth;
1524
1525 mnth = dtim->mo - 1L;
1526 mn1 = dtim->mn/10L;
1527 mn2 = dtim->mn - (mn1*10);
1528 hr1 = dtim->hr/10L;
1529 hr2 = dtim->hr - (hr1*10);
1530 dy1 = dtim->dy/10L;
1531 dy2 = dtim->dy - (dy1*10);
1532 if (tinc==1) {
1533 sprintf(ch,"%04li",dtim->yr);
1534 }
1535 else if (tinc==2) {
1536 if (dtim->yr==9999L) {
1537 sprintf(ch,"%s",monc[mnth]);
1538 } else {
1539 sprintf(ch,"%s%04li",monc[mnth],dtim->yr);
1540 }
1541 }
1542 else if (tinc==3) {
1543 sprintf(ch,"%i%i%s%04li",dy1,dy2,monc[mnth],dtim->yr);
1544 }
1545 else if (tinc==4) {
1546 sprintf(ch,"%i%iZ%i%i%s%04li",hr1,hr2,dy1,dy2,
1547 monc[mnth],dtim->yr);
1548 }
1549 else if (tinc==5) {
1550 sprintf(ch,"%i%i:%i%iZ%i%i%s%04li",hr1,hr2,mn1,mn2,dy1,dy2,
1551 monc[mnth],dtim->yr);
1552 }
1553 else sprintf(ch,"???");
1554 len=0;
1555 while (ch[len]) len++;
1556 return (len);
1557 }
1558
1559 /* Compare two strings given the length. */
1560 /* Return 0 if the string match, otherwise 1. */
1561
cmpch(char * str1,char * str2,int len)1562 int cmpch (char *str1, char *str2, int len) {
1563
1564 while (len>0) {
1565 len--;
1566 if (*(str1+len) != *(str2+len)) return (1);
1567 }
1568 return (0);
1569 }
1570
1571 /* Free anything hung off a gastat block */
1572
gafree(struct gastat * pst)1573 void gafree (struct gastat *pst) {
1574
1575 if (pst->type == 1) {
1576 gagfre (pst->result.pgr);
1577 pst->result.pgr=NULL;
1578 } else {
1579 gasfre (pst->result.stn);
1580 pst->result.stn=NULL;
1581 }
1582 }
1583
gagfre(struct gagrid * pgr)1584 void gagfre (struct gagrid *pgr) {
1585 if (pgr==NULL) return;
1586 if (pgr->alocf) {
1587 free(pgr->ivals);
1588 free(pgr->jvals);
1589 }
1590 if (pgr->idim>-1) free (pgr->grid);
1591 free (pgr);
1592 }
1593
gasfre(struct gastn * stn)1594 void gasfre (struct gastn *stn) {
1595 int i;
1596 if (stn==NULL) return;
1597 if (stn->tvals) free (stn->tvals);
1598 if (stn->rpt) {
1599 for (i=0; i<BLKNUM; i++) {
1600 if (stn->blks[i] != NULL) free (stn->blks[i]);
1601 }
1602 }
1603 free (stn);
1604 }
1605
1606
1607 /* Expand file names prefixed with '^' from data descriptor
1608 files */
1609
fnmexp(char * out,char * in1,char * in2)1610 void fnmexp (char *out, char *in1, char *in2) {
1611 char *pos, *ch, envv[20], *envr;
1612 int i,j;
1613
1614 if (*in1=='$') {
1615 in1++;
1616 i = 0;
1617 while (*in1!='/' && *in1!='\0' && i<16) {
1618 envv[i] = *in1;
1619 i++; in1++;
1620 }
1621 envv[i] = '\0';
1622 envr = gxgsym(envv);
1623 if (envr) {
1624 i = 0; j = 0;
1625 while (*(envr+j)) {
1626 *(out+i) = *(envr+j);
1627 i++; j++;
1628 }
1629 while (*in1!='\0' && *in1!=' ' && *in1!='\n') {
1630 *(out+i) = *in1;
1631 i++; in1++;
1632 }
1633 *(out+i) = '\0';
1634 }
1635 return;
1636 }
1637 ch = in2;
1638 pos=NULL;
1639 while (*ch!='\0' && *ch!=' ' && *ch!='\n') {
1640 if (*ch=='/') pos=ch;
1641 ch++;
1642 }
1643 if (pos) pos++;
1644 while (pos!=NULL && in2<pos) {
1645 *out = *in2;
1646 out++; in2++;
1647 }
1648 in1++;
1649 while (*in1!='\0' && *in1!=' ' && *in1!='\n') {
1650 *out = *in1;
1651 out++; in1++;
1652 }
1653 *out = '\0';
1654 }
1655
1656 /* Given a file name template and a dt structure, fill in to get the file name */
1657
gafndt(char * fn,struct dt * dtim,struct dt * dtimi,float * vals,struct gachsub * pch1st,int t)1658 char *gafndt (char *fn, struct dt *dtim, struct dt *dtimi, float *vals, struct gachsub *pch1st, int t) {
1659 struct gachsub *pchsub;
1660 struct dt stim;
1661 int len,olen,iv,tdif;
1662 char *fnout, *in, *out, *work, *in2, *out2;
1663
1664 olen = 0;
1665 while (*(fn+olen)) olen++;
1666 olen+=5;
1667 fnout = (char *)malloc(olen);
1668 if (fnout==NULL) return (NULL);
1669
1670 in = fn;
1671 out = fnout;
1672
1673 while (*in) {
1674 pchsub = pch1st;
1675 if (*in=='%' && *(in+1)=='i') {
1676 if (*(in+2)=='x' && *(in+3)=='1') { /* x: decades */
1677 sprintf (out,"%i",dtimi->yr/10);
1678 while (*out) out++;
1679 in+=4;
1680 } else if (*(in+2)=='x' && *(in+3)=='3') {
1681 sprintf (out,"%03i",dtimi->yr/10);
1682 out+=3; in+=4;
1683 } else if (*(in+2)=='y' && *(in+3)=='2') {
1684 iv = dtimi->yr/100;
1685 iv = dtimi->yr - iv*100;
1686 sprintf (out,"%02i",iv);
1687 out+=2; in+=4;
1688 } else if (*(in+2)=='y' && *(in+3)=='4') {
1689 sprintf (out,"%04i",dtimi->yr);
1690 out+=4; in+=4;
1691 } else if (*(in+2)=='m' && *(in+3)=='1') {
1692 sprintf (out,"%i",dtimi->mo);
1693 while (*out) out++;
1694 in+=4;
1695 } else if (*(in+2)=='m' && *(in+3)=='2') {
1696 sprintf (out,"%02i",dtimi->mo);
1697 out+=2; in+=4;
1698 } else if (*(in+2)=='m' && *(in+3)=='h') {
1699 if (dtimi->dy < 16) *out='a';
1700 else *out = 'b';
1701 out+=1; in+=4;
1702 } else if (*(in+2)=='m' && *(in+3)=='H') {
1703 if (dtimi->dy < 16) *out='A';
1704 else *out = 'B';
1705 out+=1; in+=4;
1706 } else if (*(in+2)=='m' && *(in+3)=='c') {
1707 *out = *(mons[dtimi->mo-1]);
1708 *(out+1) = *(mons[dtimi->mo-1]+1);
1709 *(out+2) = *(mons[dtimi->mo-1]+2);
1710 out+=3; in+=4;
1711 } else if (*(in+2)=='d' && *(in+3)=='1') {
1712 sprintf (out,"%i",dtimi->dy);
1713 while (*out) out++;
1714 in+=4;
1715 } else if (*(in+2)=='d' && *(in+3)=='2') {
1716 sprintf (out,"%02i",dtimi->dy);
1717 out+=2; in+=4;
1718 } else if (*(in+2)=='h' && *(in+3)=='1') {
1719 sprintf (out,"%i",dtimi->hr);
1720 while (*out) out++;
1721 in+=4;
1722 } else if (*(in+2)=='h' && *(in+3)=='2') {
1723 sprintf (out,"%02i",dtimi->hr);
1724 out+=2; in+=4;
1725 /*mf added 95082800 mf*/
1726 } else if (*(in+2)=='h' && *(in+3)=='3') {
1727 sprintf (out,"%03i",dtimi->hr);
1728 out+=3; in+=4;
1729 } else if (*(in+2)=='n' && *(in+3)=='2') {
1730 sprintf (out,"%02i",dtimi->mn);
1731 out+=2; in+=4;
1732 } else {
1733 *out = *in;
1734 in++; out++;
1735 }
1736 } else if (*in=='%' && *(in+1)=='x' && *(in+2)=='1') { /* x: decades */
1737 sprintf (out,"%i",dtim->yr/10);
1738 while (*out) out++;
1739 in+=3;
1740 } else if (*in=='%' && *(in+1)=='x' && *(in+2)=='3') {
1741 sprintf (out,"%03i",dtim->yr/10);
1742 out+=3; in+=3;
1743 } else if (*in=='%' && *(in+1)=='y' && *(in+2)=='2') {
1744 iv = dtim->yr/100;
1745 iv = dtim->yr - iv*100;
1746 sprintf (out,"%02i",iv);
1747 out+=2; in+=3;
1748 } else if (*in=='%' && *(in+1)=='y' && *(in+2)=='4') {
1749 sprintf (out,"%04i",dtim->yr);
1750 out+=4; in+=3;
1751 } else if (*in=='%' && *(in+1)=='m' && *(in+2)=='1') {
1752 sprintf (out,"%i",dtim->mo);
1753 while (*out) out++;
1754 in+=3;
1755 } else if (*in=='%' && *(in+1)=='m' && *(in+2)=='2') {
1756 sprintf (out,"%02i",dtim->mo);
1757 out+=2; in+=3;
1758 } else if (*in=='%' && *(in+1)=='m' && *(in+2)=='h') {
1759 if (dtim->dy < 16) *out='a';
1760 else *out = 'b';
1761 out+=1; in+=3;
1762 } else if (*in=='%' && *(in+1)=='m' && *(in+2)=='H') {
1763 if (dtim->dy < 16) *out='A';
1764 else *out = 'B';
1765 out+=1; in+=3;
1766 } else if (*in=='%' && *(in+1)=='m' && *(in+2)=='c') {
1767 *out = *(mons[dtim->mo-1]);
1768 *(out+1) = *(mons[dtim->mo-1]+1);
1769 *(out+2) = *(mons[dtim->mo-1]+2);
1770 out+=3; in+=3;
1771 } else if (*in=='%' && *(in+1)=='d' && *(in+2)=='1') {
1772 sprintf (out,"%i",dtim->dy);
1773 while (*out) out++;
1774 in+=3;
1775 } else if (*in=='%' && *(in+1)=='d' && *(in+2)=='2') {
1776 sprintf (out,"%02i",dtim->dy);
1777 out+=2; in+=3;
1778 } else if (*in=='%' && *(in+1)=='h' && *(in+2)=='1') {
1779 sprintf (out,"%i",dtim->hr);
1780 while (*out) out++;
1781 in+=3;
1782 } else if (*in=='%' && *(in+1)=='h' && *(in+2)=='2') {
1783 sprintf (out,"%02i",dtim->hr);
1784 out+=2; in+=3;
1785 /*mf added 95082800 mf*/
1786 } else if (*in=='%' && *(in+1)=='h' && *(in+2)=='3') {
1787 sprintf (out,"%03i",dtim->hr);
1788 out+=3; in+=3;
1789 } else if (*in=='%' && *(in+1)=='n' && *(in+2)=='2') {
1790 sprintf (out,"%02i",dtim->mn);
1791 out+=2; in+=3;
1792 } else if (*in=='%' && *(in+1)=='f' && *(in+2)=='2') {
1793 stim.yr = (int)(*vals+0.1);
1794 stim.mo = (int)(*(vals+1)+0.1);
1795 stim.dy = (int)(*(vals+2)+0.1);
1796 stim.hr = (int)(*(vals+3)+0.1);
1797 stim.mn = (int)(*(vals+4)+0.1);
1798 tdif = timdif(&stim,dtim);
1799 tdif = (tdif+30)/60;
1800 if (tdif<99) sprintf (out,"%02i",tdif);
1801 else sprintf (out,"%i",tdif);
1802 while (*out) out++;
1803 in+=3;
1804 /*mf added 95082800 mf*/
1805 } else if (*in=='%' && *(in+1)=='f' && *(in+2)=='3') {
1806 stim.yr = (int)(*vals+0.1);
1807 stim.mo = (int)(*(vals+1)+0.1);
1808 stim.dy = (int)(*(vals+2)+0.1);
1809 stim.hr = (int)(*(vals+3)+0.1);
1810 stim.mn = (int)(*(vals+4)+0.1);
1811 tdif = timdif(&stim,dtim);
1812 tdif = (tdif+30)/60;
1813 if (tdif<999) sprintf (out,"%03i",tdif);
1814 else sprintf (out,"%i",tdif);
1815 while (*out) out++;
1816 in+=3;
1817 } else if (*in=='%' && *(in+1)=='c' && *(in+2)=='h') {
1818 while (pchsub) {
1819 if (t>=pchsub->t1 && (pchsub->t2 == -99 || t<=pchsub->t2) ) {
1820 len = wrdlen(pchsub->ch); /* Reallocate output string */
1821 olen += len;
1822 work = (char *)malloc(olen);
1823 if (work==NULL) {
1824 free (fnout);
1825 return (NULL);
1826 }
1827 in2 = fnout;
1828 out2 = work;
1829 while (in2!=out) {
1830 *out2 = *in2;
1831 in2++; out2++;
1832 }
1833 free (fnout);
1834 fnout = work;
1835 out = out2;
1836 getwrd(out,pchsub->ch,len);
1837 out += len;
1838 break;
1839 }
1840 pchsub = pchsub->forw;
1841 }
1842 in+=3;
1843 } else {
1844 *out = *in;
1845 in++; out++;
1846 }
1847 }
1848 *out = '\0';
1849 return (fnout);
1850 }
1851
1852 /* Byte swap requested number of 4 byte elements */
1853
gabswp(void * r,int cnt)1854 void gabswp (void *r, int cnt) {
1855 int i;
1856 char *ch1,*ch2,*ch3,*ch4,cc1,cc2;
1857
1858 ch1 = (char *)r;
1859 ch2 = ch1+1;
1860 ch3 = ch2+1;
1861 ch4 = ch3+1;
1862 for (i=0; i<cnt; i++) {
1863 cc1 = *ch1;
1864 cc2 = *ch2;
1865 *ch1 = *ch4;
1866 *ch2 = *ch3;
1867 *ch3 = cc2;
1868 *ch4 = cc1;
1869 ch1+=4; ch2+=4; ch3+=4; ch4+=4;
1870 }
1871 }
1872
1873 /* joew 071902 Byte swap a single element with the width given */
ganbswp(char * buf,int cnt)1874 void ganbswp(char *buf, int cnt) {
1875 int i, j;
1876 char tmp;
1877 for (i = 0, j = cnt - 1; i < cnt / 2; i++, j--) {
1878 tmp = buf[i];
1879 buf[i] = buf[j];
1880 buf[j] = tmp;
1881 }
1882 }
1883
1884 /* Byte swap a report header from a station data file */
1885
gahswp(struct rpthdr * hdr)1886 void gahswp (struct rpthdr *hdr) {
1887 gabswp((float *)(&(hdr->lat)),5);
1888 }
1889
1890 /* Return day of week for date/time 0=sunday, 6=saturday */
1891
dayweek(struct dt * dtime)1892 int dayweek (struct dt *dtime) {
1893 struct dt anch;
1894 int i,j;
1895 if (dtime->yr<1950 || dtime->yr>2020) return(7);
1896 anch.yr = 1950;
1897 anch.mo = 1;
1898 anch.dy = 1;
1899 anch.hr = 0;
1900 anch.mn = 0;
1901 i = timdif(&anch,dtime);
1902 i = i/1440;
1903 j = i/7;
1904 i = i - j*7;
1905 return(i);
1906 }
1907
1908 /*
1909 * convert an IMB float to single precision number v1.0
1910 *
1911 * Wesley Ebisuzaki
1912 */
1913
ibm2flt(unsigned char * ibm)1914 float ibm2flt(unsigned char *ibm) {
1915
1916 int positive, power;
1917 unsigned int abspower;
1918 long int mant;
1919 double value, exp;
1920
1921 positive = (ibm[0] & 0x80) == 0;
1922 mant = (ibm[1] << 16) + (ibm[2] << 8) + ibm[3];
1923 power = (int) (ibm[0] & 0x7f) - 64;
1924 abspower = power > 0 ? power : -power;
1925
1926
1927 /* calc exp */
1928 exp = 16.0;
1929 value = 1.0;
1930 while (abspower) {
1931 if (abspower & 1) {
1932 value *= exp;
1933 }
1934 exp = exp * exp;
1935 abspower >>= 1;
1936 }
1937
1938 if (power < 0) value = 1.0 / value;
1939 value = value * mant / 16777216.0;
1940 if (positive == 0) value = -value;
1941 return (float)value;
1942 }
1943
1944 /*
1945 * convert a float to an IBM single precision number v1.0
1946 *
1947 * Wesley Ebisuzaki
1948 *
1949 * doesn't handle subnormal numbers
1950 */
1951
flt2ibm(float x,unsigned char * ibm)1952 int flt2ibm(float x, unsigned char *ibm) {
1953
1954 int sign, exp, i;
1955 double mant;
1956
1957 if (x == 0.0) {
1958 ibm[0] = ibm[1] = ibm[2] = ibm[3] = 0;
1959 return 0;
1960 }
1961
1962 /* sign bit */
1963 if (x < 0.0) {
1964 sign = 128;
1965 x = -x;
1966 }
1967 else sign = 0;
1968
1969 mant = frexp((double) x, &exp);
1970
1971 /* round up by adding 2**-24 */
1972 /* mant = mant + 1.0/16777216.0; */
1973
1974 if (mant >= 1.0) {
1975 mant = 0.5;
1976 exp++;
1977 }
1978 while (exp & 3) {
1979 mant *= 0.5;
1980 exp++;
1981 }
1982
1983 exp = exp/4 + 64;
1984
1985 if (exp < 0) {
1986 fprintf(stderr,"underflow in flt2ibm\n");
1987 ibm[0] = ibm[1] = ibm[2] = ibm[3] = 0;
1988 return 0;
1989 }
1990 if (exp > 127) {
1991 fprintf(stderr,"overflow in flt2ibm\n");
1992 ibm[0] = sign | 127;
1993 ibm[1] = ibm[2] = ibm[3] = 255;
1994 return -1;
1995 }
1996
1997 /* normal number */
1998
1999 ibm[0] = sign | exp;
2000
2001 mant = mant * 256.0;
2002 i = floor(mant);
2003 mant = mant - i;
2004 ibm[1] = i;
2005
2006 mant = mant * 256.0;
2007 i = floor(mant);
2008 mant = mant - i;
2009 ibm[2] = i;
2010
2011 ibm[3] = floor(mant*256.0);
2012
2013 return 0;
2014 }
2015
2016 /* wesley ebisuzaki v0.1
2017 *
2018 * takes 4 byte character string (single precision ieee big-endian)
2019 * and returns a float
2020 *
2021 * doesn't handle NAN, infinity and any other funny stuff in ieee
2022 *
2023 * ansi C
2024 */
2025
ieee2flt(unsigned char * ieee)2026 float ieee2flt(unsigned char *ieee) {
2027 double fmant;
2028 int exp;
2029
2030 if (ieee[0] == 0 && ieee[1] == 0 && ieee[2] == 0 && ieee[3] == 0)
2031 return (float) 0.0;
2032
2033 exp = ((ieee[0] & 127) << 1) + (ieee[1] >> 7);
2034 fmant = (double) ((int) ieee[3] + (int) (ieee[2] << 8) +
2035 (int) ((ieee[1] | 128) << 16));
2036 if (ieee[0] & 128) fmant = -fmant;
2037 return (float) (ldexp(fmant, (int) (exp - 128 - 22)));
2038 }
2039
2040
2041 /*
2042 * convert a float to an ieee single precision number v1.1
2043 * (big endian)
2044 * Wesley Ebisuzaki
2045 *
2046 * bugs: doesn't handle subnormal numbers
2047 * bugs: assumes length of integer >= 25 bits
2048 */
2049
flt2ieee(float x,unsigned char * ieee)2050 int flt2ieee(float x, unsigned char *ieee) {
2051
2052 int sign, exp;
2053 unsigned int umant;
2054 double mant;
2055
2056 if (x == 0.0) {
2057 ieee[0] = ieee[1] = ieee[2] = ieee[3] = 0;
2058 return 0;
2059 }
2060
2061 /* sign bit */
2062 if (x < 0.0) {
2063 sign = 128;
2064 x = -x;
2065 }
2066 else sign = 0;
2067 mant = frexp((double) x, &exp);
2068
2069 /* 2^24 = 16777216 */
2070
2071 umant = mant * 16777216 + 0.5;
2072 if (umant >= 16777216) {
2073 umant = umant / 2;
2074 exp++;
2075 }
2076 /* bit 24 should be a 1 .. not used in ieee format */
2077
2078 exp = exp - 1 + 127;
2079
2080 if (exp < 0) {
2081 /* signed zero */
2082 ieee[0] = sign;
2083 ieee[1] = ieee[2] = ieee[3] = 0;
2084 return 0;
2085 }
2086 if (exp > 255) {
2087 /* signed infinity */
2088 ieee[0] = sign + 127;
2089 ieee[1] = 128;
2090 ieee[2] = ieee[3] = 0;
2091 return 0;
2092 }
2093 /* normal number */
2094
2095 ieee[0] = sign + (exp >> 1);
2096
2097 ieee[3] = umant & 255;
2098 ieee[2] = (umant >> 8) & 255;
2099 ieee[1] = ((exp & 1) << 7) + ((umant >> 16) & 127);
2100 return 0;
2101 }
2102
be_int2int(unsigned char * be_int)2103 int be_int2int(unsigned char *be_int) {
2104
2105 int sign;
2106 unsigned long n;
2107
2108 sign = (be_int[0] & 128);
2109 n = (unsigned int) ((be_int[0] & 127) << 24) +
2110 (unsigned int) (be_int[1] << 16) +
2111 (unsigned int) (be_int[2] << 8) +
2112 (unsigned int) be_int[3];
2113
2114 /* positive number */
2115
2116 if (sign == 0) {
2117 return (n <= INT_MAX) ? (int) n : (int) INT_MAX;
2118 }
2119
2120 /* negative number */
2121 n = (n ^ 0x7fffffff);
2122
2123 if (n <= (unsigned) -(INT_MIN+1)) {
2124 return -1 - (int) n;
2125 }
2126 return INT_MIN;
2127 }
2128
2129 /* Copies indicated scaling info into newly allocated
2130 floating point array. args:
2131
2132 vals -- input scaling array
2133 lin -- input is linear or levels
2134 dir -- direction of scaling info. 0 for gr to ab;
2135 1 for ab to gr
2136 dim -- dimension the scaling info is for
2137
2138 lin, dir, and dim are provided solely to figure out how
2139 many values are to be copied. This assumes knowledge
2140 of how the various scaling items are set up. */
2141
2142
cpscal(float * vals,int lin,int dir,int dim)2143 float *cpscal (float *vals, int lin, int dir, int dim) {
2144 int i,num;
2145 float *vvv;
2146
2147 if (dim<0) return (NULL);
2148 if (dim==3) num = 8;
2149 else {
2150 if (lin==1) num = 3;
2151 else num = (int)(*vals+0.5) + 5;
2152 }
2153 vvv = (float *)malloc(sizeof(float)*num);
2154 if (vvv==NULL) return (NULL);
2155 for (i=0; i<num; i++) *(vvv+i) = *(vals+i);
2156 return (vvv);
2157 }
2158
2159 /* handle var name of the form longnm=>abbrv
2160 or just the abbrv with no long name */
2161
getvnm(struct gavar * pvar,char * mrec)2162 int getvnm (struct gavar *pvar, char *mrec) {
2163 int ib,i,j,k,len,flag;
2164 char *ch;
2165
2166 ib = 0;
2167 while (*(mrec+ib)==' ') ib++;
2168
2169 if (*(mrec+ib)=='\0' || *(mrec+ib)=='\n') return(1);
2170
2171 /* Scan for the '=>' string */
2172 len = 0;
2173 i = ib;
2174 flag = 0;
2175
2176 while (1) {
2177 if (*(mrec+i)==' ' || *(mrec+i)=='\0' || *(mrec+i)=='\n') break;
2178 if (*(mrec+i)=='=' && *(mrec+i+1)=='>') {
2179 flag = 1;
2180 break;
2181 }
2182 len++ ; i++;
2183 }
2184
2185 if (flag) {
2186 ch = (char *)malloc(len+1);
2187 if (ch==NULL) return(1);
2188 for (j=ib; j<i; j++) {
2189 k = j-ib;
2190 *(ch+k) = *(mrec+j);
2191 }
2192 *(ch+len) = '\0';
2193 pvar->longnm = ch;
2194 i+=2;
2195 } else {
2196 i = 0;
2197 }
2198
2199 if (*(mrec+i)=='\n' || *(mrec+i)=='\0') return (1);
2200
2201 getwrd (pvar->abbrv, mrec+i, 15);
2202 lowcas(pvar->abbrv);
2203
2204 /* Check if 1st character is lower-case alphabetic */
2205 if (islower(*(pvar->abbrv))) return(0);
2206 else return (1);
2207 }
2208
2209 #ifndef HAVE_FSEEKO
2210
fseeko(FILE * stream,off_t offset,int whence)2211 int fseeko(FILE *stream, off_t offset, int whence) {
2212 fseek(stream, (long)offset, whence);
2213 }
2214
ftello(FILE * stream)2215 off_t ftello(FILE *stream) {
2216 return (off_t)ftell(stream);
2217 }
2218 #endif
2219
2220
2221
2222