3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996 Dennis F. Flanigan Jr. of Applied Research Corporation,
5 Landover, Maryland, a NASA/GSFC on-site contractor.
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Library General Public
9 License as published by the Free Software Foundation; either
10 version 2 of the License, or (at your option) any later version.
12 This library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 Library General Public License for more details.
17 You should have received a copy of the GNU Library General Public
18 License along with this library; if not, write to the Free
19 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22 * Interpolation functions
24 * Dennis Flanigan, Jr.
27 * Finished testing interpolation routine.
30 * Finished up bilinear interpolation routine.
31 * Routine calculates values using four surrounding
32 * values with the same slant range.
35 * Added bilinear interpolation value routine. Rewrote
36 * get_surrounding sweeps so that it returns Sweeps instead of
40 * Added internal ray searching routine. Replaces
41 * RSL_get_next_closest_ray.
44 * Added internal sweep searching routine designed for use
45 * by interpolation routines. Replaced RSL_get_next_closest_sweep.
48 * This file was created. Started adding internal routines
49 * needed for calculating distences between two points
66 extern double hypot(double, double);
68 /***************************************************/
73 /***************************************************/
74 void get_xyz_coord(double *x,double *y,double *z,
75 double range,double azim,double elev)
77 /* Return x,y,z coordinates given range, azimuth angle and
78 * elevation angle. Memory allocation for x,y and z are
79 * not provided by this routine!
81 double azim_rad,elev_rad;
83 # define M_PI 3.14159265358979323846
84 azim_rad = azim * M_PI / 180.0;
85 elev_rad = elev * M_PI / 180.0;
87 *x = cos(elev_rad) * cos(azim_rad) * range;
88 *y = cos(elev_rad) * sin(azim_rad) * range;
89 *z = sin(elev_rad) * range;
92 /***************************************************/
97 /***************************************************/
98 float get_dist(float r1,float a1,float e1,
99 float r2,float a2,float e2)
101 /* Give two points described by range, azimuth angle
102 * and elevation angle, return the distence between
105 double x1,y1,z1,x2,y2,z2;
107 get_xyz_coord(&x1,&y1,&z1,(double)r1,(double)a1,(double)e1);
108 get_xyz_coord(&x2,&y2,&z2,(double)r2,(double)a2,(double)e2);
110 return (float) sqrt(pow(x1 - x2,2) + pow(y1 - y2,2) + pow(z1 - z2,2));
113 /***************************************************/
115 /* get_surrounding_sweeps */
118 /***************************************************/
119 void get_surrounding_sweep(Sweep **below,Sweep **above, Volume *v,
122 /* Return the pointers of the sweeps that are above and
123 * and below the elevation angle in the parameter list.
125 * Assume at least one non-NULL sweep exist in Volume.
127 * A value of NULL is set to above or below in cases
128 * where there is no sweep above or below.
132 /* look for above index first */
134 while(a < v->h.nsweeps)
136 if(v->sweep[a] != NULL)
138 if(v->sweep[a]->h.elev >= elev)
140 *above = v->sweep[a];
147 /* Was above ever set ?
148 * If not, set above to counter.
150 if(a == v->h.nsweeps)
155 /* Look for index just below elev */
159 if(v->sweep[a] != NULL)
161 if(v->sweep[a]->h.elev <= elev)
163 *below = v->sweep[a];
170 /* Was a below found ? */
179 /******************************************
183 * Dennis Flanigan,Jr. 4/29/95 *
184 ******************************************/
185 double dir_angle_diff(float x,float y)
187 /* returns difference between angles x and y. Returns
188 * positive value if y > x, negitive value if
195 while (d >= 180) d = -1 * (360 - d);
196 while (d < -180) d = (360 + d);
201 /***************************************************/
203 /* get_surrounding_rays */
206 /***************************************************/
207 void get_surrounding_ray(Ray **ccwise,Ray **cwise,Sweep *s,
210 /* Return the pointers to the rays that are counterclockwise
211 * and clockwise to the ray_angle in the parameter list.
212 * Memory space for the variables ccwise and cwise
213 * is not allocated by this routine.
215 * Assume at least two rays exist and that a hash
216 * table has been created.
218 * Will need to add test for first ray and last
219 * ray if this routine is going to be used for
224 Hash_table *hash_table;
226 Azimuth_hash *closest;
228 /* Find hash index close to hindex we want. This will
229 * used as a starting point for a search to the closest
232 hash_table = hash_table_for_sweep(s);
233 if (hash_table == NULL) return; /* Nada. */
235 hindex = hash_bin(hash_table,ray_angle);
237 /* Find hash entry with closest Ray */
238 closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
240 close_diff = dir_angle_diff(ray_angle,closest->ray->h.azimuth);
244 /* Closest ray is counterclockwise to ray_angle */
245 *ccwise = closest->ray;
246 *cwise = closest->ray_high->ray;
250 /* Closest ray is clockwise to ray_angle. */
251 *cwise = closest->ray;
252 *ccwise = closest->ray_low->ray;
257 /***************************************************/
261 /***************************************************/
262 double from_dB(double db)
264 return pow(10,db/10.0);
267 double to_dB(double value)
269 return 10.0 * log10(value);
272 /***************************************************/
274 /* get_linear_value_from_sweep */
276 /***************************************************/
277 double get_linear_value_from_sweep(Sweep *sweep,float srange,float azim,
280 float ccw_db_value,cw_db_value;
281 double ccw_value,cw_value,value;
283 Ray *ccw_ray,*cw_ray;
285 get_surrounding_ray(&ccw_ray,&cw_ray,sweep,azim);
287 /* Assume that ccw_ray and cw_ray will be non_NULL */
289 if((azim - ccw_ray->h.azimuth) > limit)
295 ccw_db_value = RSL_get_value_from_ray(ccw_ray,srange);
297 /* Check to make sure this is a valid value */
298 if (ccw_db_value == BADVAL)
304 ccw_value = from_dB(ccw_db_value);
308 if((cw_ray->h.azimuth - azim) > limit)
314 cw_db_value = RSL_get_value_from_ray(cw_ray,srange);
315 /* Check to make sure this is a valid value */
316 if (cw_db_value == BADVAL)
322 cw_value = from_dB(cw_db_value);
326 if((cw_value != -1) && (ccw_value != -1))
328 /* Both the clockwise ray and the counterclockwise
331 delta_angle = angle_diff(ccw_ray->h.azimuth,cw_ray->h.azimuth);
333 value=((angle_diff(azim,cw_ray->h.azimuth)/delta_angle)*ccw_value)
334 + ((angle_diff(azim,ccw_ray->h.azimuth)/delta_angle) * cw_value);
336 else if((cw_value == -1) && (ccw_value == -1))
338 /* Neither ray is valid. */
341 else if(cw_value != -1)
343 /* cw_ray is only ray that is within limit. */
348 /* ccw_ray is only ray that is within limit. */
355 /***************************************************/
357 /* RSL_get_linear_value */
360 /***************************************************/
361 float RSL_get_linear_value(Volume *v,float srange,float azim,
362 float elev,float limit)
364 /* Compute bilinear value from four surrounding values
365 * in Volume v. The four surrounding values
366 * are at a constant range.
368 * Limit is an angle used only to reject values
369 * in the azimuth plane.
373 double up_value, down_value;
375 Sweep *up_sweep,*down_sweep;
377 get_surrounding_sweep(&down_sweep,&up_sweep,v,elev);
379 /* Calculate interpolated value in sweep above
388 up_value = get_linear_value_from_sweep(up_sweep,srange,azim,limit);
391 /* Calculate interpolated value in sweep below requested point.
393 if(down_sweep == NULL)
399 down_value = get_linear_value_from_sweep(down_sweep,srange,azim,limit);
401 /* Using the interpolated values calculated at the elevation
402 * angles in the above and below sweeps, interpolate a value
403 * for the elvation angle at the requested point.
405 if((up_value != -1) && (down_value != -1))
407 delta_angle = angle_diff(up_sweep->h.elev,down_sweep->h.elev);
409 value =((angle_diff(elev,up_sweep->h.elev)/delta_angle) * down_value) +
410 ((angle_diff(elev,down_sweep->h.elev)/delta_angle) * up_value);
412 else if((up_value == -1) && (down_value == -1))
416 else if(up_value != -1)
425 /* Convert back to dB value and return. */
428 db_value = (float)to_dB(value);