1 /*
2 * (C) Copyright 2005- ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10 /**************************************
11 * Enrico Fucile
12 **************************************/
13
14
15 #include "grib_api_internal.h"
16 #include <math.h>
17
18 /*
19 This is used by make_class.pl
20
21 START_CLASS_DEF
22 CLASS = iterator
23 SUPER = grib_iterator_class_regular
24 IMPLEMENTS = init
25 END_CLASS_DEF
26
27 */
28
29 /* START_CLASS_IMP */
30
31 /*
32
33 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
34 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
35 or edit "iterator.class" and rerun ./make_class.pl
36
37 */
38
39
40 static void init_class(grib_iterator_class*);
41
42 static int init(grib_iterator* i, grib_handle*, grib_arguments*);
43
44
45 typedef struct grib_iterator_gaussian
46 {
47 grib_iterator it;
48 /* Members defined in gen */
49 long carg;
50 const char* missingValue;
51 /* Members defined in regular */
52 double* las;
53 double* los;
54 long Ni;
55 long Nj;
56 long iScansNegatively;
57 long isRotated;
58 double angleOfRotation;
59 double southPoleLat;
60 double southPoleLon;
61 long jPointsAreConsecutive;
62 long disableUnrotate;
63 /* Members defined in gaussian */
64 } grib_iterator_gaussian;
65
66 extern grib_iterator_class* grib_iterator_class_regular;
67
68 static grib_iterator_class _grib_iterator_class_gaussian = {
69 &grib_iterator_class_regular, /* super */
70 "gaussian", /* name */
71 sizeof(grib_iterator_gaussian), /* size of instance */
72 0, /* inited */
73 &init_class, /* init_class */
74 &init, /* constructor */
75 0, /* destructor */
76 0, /* Next Value */
77 0, /* Previous Value */
78 0, /* Reset the counter */
79 0, /* has next values */
80 };
81
82 grib_iterator_class* grib_iterator_class_gaussian = &_grib_iterator_class_gaussian;
83
84
init_class(grib_iterator_class * c)85 static void init_class(grib_iterator_class* c)
86 {
87 c->next = (*(c->super))->next;
88 c->previous = (*(c->super))->previous;
89 c->reset = (*(c->super))->reset;
90 c->has_next = (*(c->super))->has_next;
91 }
92 /* END_CLASS_IMP */
93
94 static void binary_search(double xx[], const unsigned long n, double x, unsigned long* j);
95
init(grib_iterator * i,grib_handle * h,grib_arguments * args)96 static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
97 {
98 grib_iterator_gaussian* self = (grib_iterator_gaussian*)i;
99
100 double* lats;
101 double laf; /* latitude of first point in degrees */
102 double lal; /* latitude of last point in degrees */
103 long trunc; /* number of parallels between a pole and the equator */
104 long lai;
105 long jScansPositively = 0;
106 int size;
107 double start;
108 unsigned long istart = 0;
109 int ret = GRIB_SUCCESS;
110
111 const char* latofirst = grib_arguments_get_name(h, args, self->carg++);
112 const char* latoflast = grib_arguments_get_name(h, args, self->carg++);
113 const char* numtrunc = grib_arguments_get_name(h, args, self->carg++);
114 const char* s_jScansPositively = grib_arguments_get_name(h, args, self->carg++);
115
116 if ((ret = grib_get_double_internal(h, latofirst, &laf)))
117 return ret;
118 if ((ret = grib_get_double_internal(h, latoflast, &lal)))
119 return ret;
120 if ((ret = grib_get_long_internal(h, numtrunc, &trunc)))
121 return ret;
122 if ((ret = grib_get_long_internal(h, s_jScansPositively, &jScansPositively)))
123 return ret;
124
125 start = laf;
126
127 size = trunc * 2;
128
129 lats = (double*)grib_context_malloc(h->context, size * sizeof(double));
130
131 ret = grib_get_gaussian_latitudes(trunc, lats);
132
133 if (ret != GRIB_SUCCESS) {
134 grib_context_log(h->context, GRIB_LOG_ERROR, "error %d calculating gaussian points", ret);
135 return ret;
136 }
137 /*
138 for(loi=(trunc*2)-1;loi>=0;loi--)
139 if(fabs(lats[loi] - lal) < glatPrecision) break;
140 for(j=(trunc*2)-1;j>0;j--) {
141 if(fabs(lats[j] - laf) < glatPrecision) break;
142 }
143 */
144
145 binary_search(lats, size - 1, start, &istart);
146 Assert(istart < size);
147
148 if (jScansPositively) {
149 for (lai = 0; lai < self->Nj; lai++) {
150 self->las[lai] = lats[istart--];
151 /*if (istart<0) istart=size-1; this condition is always FALSE -- 'istart' is unsigned long */
152 }
153 }
154 else {
155 for (lai = 0; lai < self->Nj; lai++) {
156 self->las[lai] = lats[istart++];
157 if (istart > size - 1)
158 istart = 0;
159 }
160 }
161
162 grib_context_free(h->context, lats);
163
164 return ret;
165 }
166
binary_search(double xx[],const unsigned long n,double x,unsigned long * j)167 static void binary_search(double xx[], const unsigned long n, double x, unsigned long* j)
168 {
169 /*This routine works only on descending ordered arrays*/
170 #define EPSILON 1e-3
171
172 unsigned long ju, jm, jl;
173 jl = 0;
174 ju = n;
175 while (ju - jl > 1) {
176 jm = (ju + jl) >> 1;
177 if (fabs(x - xx[jm]) < EPSILON) {
178 /* found something close enough. We're done */
179 *j = jm;
180 return;
181 }
182 if (x < xx[jm])
183 jl = jm;
184 else
185 ju = jm;
186 }
187 *j = jl;
188 }
189