Pileus Git - ~andy/rsl/blob - interp.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   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.
   6
   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.
  11
  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.
  16
  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.
  20 */
  21 /*
  22  * Interpolation functions
  23  *
  24  * Dennis Flanigan, Jr.
  25  *
  26  * 7/7/95
  27  * Finished testing interpolation routine.
  28  *
  29  * 7/6/95
  30  * Finished up bilinear interpolation routine.
  31  * Routine calculates values using four surrounding
  32  * values with the same slant range.
  33  *
  34  * 6/29/95
  35  * Added bilinear interpolation value routine.   Rewrote
  36  * get_surrounding sweeps so that it returns Sweeps instead of
  37  * Sweep indexes.
  38  *
  39  * 6/28/95
  40  * Added internal ray searching routine. Replaces
  41  * RSL_get_next_closest_ray.
  42  *
  43  * 6/25/95
  44  * Added internal sweep searching routine designed for use
  45  * by interpolation routines.  Replaced RSL_get_next_closest_sweep.
  46  *
  47  * 6/23/95
  48  * This file was created.  Started adding internal routines
  49  * needed for calculating distences between two points
  50  * in space.
  51  */
  52
  53 #include <stdio.h>
  54 #include <math.h>
  55 #include "rsl.h"
  56
  57
  58 #ifndef FALSE
  59 #define FALSE 0
  60 #endif
  61
  62 #ifndef TRUE
  63 #define TRUE 1
  64 #endif
  65
  66 extern double hypot(double, double);
  67
  68 /***************************************************/
  69 /*                                                 */
  70 /*          get_xyz_coord                          */
  71 /*                                                 */
  72 /*                                                 */
  73 /***************************************************/
  74 void get_xyz_coord(double *x,double *y,double *z,
  75                    double range,double azim,double elev)
  76    {
  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!
  80         */
  81    double azim_rad,elev_rad;
  82
  83 #  define M_PI          3.14159265358979323846
  84    azim_rad = azim *  M_PI / 180.0;
  85    elev_rad = elev *  M_PI / 180.0;
  86
  87    *x = cos(elev_rad) * cos(azim_rad) * range;
  88    *y = cos(elev_rad) * sin(azim_rad) * range;
  89    *z = sin(elev_rad) * range;
  90    }
  91
  92 /***************************************************/
  93 /*                                                 */
  94 /*          get_dist                               */
  95 /*                                                 */
  96 /*                                                 */
  97 /***************************************************/
  98 float get_dist(float r1,float a1,float e1,
  99                float r2,float a2,float e2)
 100    {
 101    /* Give two points described by range, azimuth angle
 102     * and elevation angle, return the distence between
 103     * the two.
 104         */
 105    double x1,y1,z1,x2,y2,z2;
 106
 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);
 109
 110    return (float) sqrt(pow(x1 - x2,2) + pow(y1 - y2,2) + pow(z1 - z2,2));
 111    }
 112
 113 /***************************************************/
 114 /*                                                 */
 115 /*          get_surrounding_sweeps                 */
 116 /*                                                 */
 117 /*                                                 */
 118 /***************************************************/
 119 void get_surrounding_sweep(Sweep **below,Sweep **above, Volume *v,
 120                            float elev)
 121    {
 122    /* Return the pointers of the sweeps that are above and
 123     * and below the elevation angle in the parameter list.
 124     *
 125     * Assume at least one non-NULL sweep exist in Volume.
 126     *
 127     * A value of NULL is set to  above or below in cases
 128     * where there is no sweep above or below.
 129     */
 130    int a;
 131
 132    /* look for above index first */
 133    a = 0;
 134    while(a < v->h.nsweeps)
 135       {
 136       if(v->sweep[a] != NULL)
 137          {
 138          if(v->sweep[a]->h.elev >= elev)
 139             {
 140             *above = v->sweep[a];
 141             break;
 142             }
 143          }
 144       a++;
 145       }
 146
 147    /* Was above ever set ?
 148     * If not, set above to counter.
 149     */
 150    if(a == v->h.nsweeps)
 151       {
 152       *above = NULL;
 153       }
 154
 155    /* Look for index just below elev */
 156    a--;
 157    while(a >= 0)
 158       {
 159       if(v->sweep[a] != NULL)
 160          {
 161          if(v->sweep[a]->h.elev <= elev)
 162             {
 163             *below = v->sweep[a];
 164             break;
 165             }
 166          }
 167       a--;
 168       }
 169
 170    /* Was a below found ? */
 171    if (a == -1)
 172       {
 173       *below = NULL;
 174       }
 175
 176    /* Bye */
 177    }
 178
 179 /******************************************
 180  *                                        *
 181  * dir_angle_diff                         *
 182  *                                        *
 183  * Dennis Flanigan,Jr. 4/29/95            *
 184  ******************************************/
 185 double dir_angle_diff(float x,float y)
 186    {
 187    /* returns difference between angles x and y.  Returns
 188     * positive value if y > x, negitive value if
 189     * y < x.
 190     *
 191     */
 192    double d;
 193    d = (double)(y - x);
 194
 195    while (d >= 180) d = -1 * (360 - d);
 196    while (d < -180) d =  (360 + d);
 197
 198    return d;
 199    }
 200
 201 /***************************************************/
 202 /*                                                 */
 203 /*          get_surrounding_rays                   */
 204 /*                                                 */
 205 /*                                                 */
 206 /***************************************************/
 207 void get_surrounding_ray(Ray **ccwise,Ray **cwise,Sweep *s,
 208                          float ray_angle)
 209    {
 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.
 214     *
 215     *  Assume at least two rays exist and that a hash
 216     *  table has been created.
 217     *
 218     *  Will need to add test for first ray and last
 219     *  ray if this routine is going to be used for
 220     *  RHI scans.
 221         */
 222    int hindex;
 223    double close_diff;
 224    Hash_table *hash_table;
 225
 226    Azimuth_hash *closest;
 227
 228    /* Find hash index close to hindex we want.  This will
 229     * used as a starting point for a search to the closest
 230     * ray.
 231     */
 232    hash_table = hash_table_for_sweep(s);
 233    if (hash_table == NULL) return; /* Nada. */
 234
 235    hindex = hash_bin(hash_table,ray_angle);
 236
 237    /* Find hash entry with closest Ray */
 238    closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
 239
 240    close_diff = dir_angle_diff(ray_angle,closest->ray->h.azimuth);
 241
 242    if(close_diff < 0)
 243       {
 244       /* Closest ray is counterclockwise to ray_angle */
 245       *ccwise = closest->ray;
 246       *cwise = closest->ray_high->ray;
 247       }
 248    else
 249       {
 250       /* Closest ray is clockwise to ray_angle. */
 251       *cwise = closest->ray;
 252       *ccwise = closest->ray_low->ray;
 253       }
 254
 255    }
 256
 257 /***************************************************/
 258 /*                                                 */
 259 /*          from_dB, to_dB                         */
 260 /*                                                 */
 261 /***************************************************/
 262 double from_dB(double db)
 263    {
 264    return pow(10,db/10.0);
 265    }
 266
 267 double to_dB(double value)
 268    {
 269    return 10.0 * log10(value);
 270    }
 271
 272 /***************************************************/
 273 /*                                                 */
 274 /*        get_linear_value_from_sweep              */
 275 /*                                                 */
 276 /***************************************************/
 277 double get_linear_value_from_sweep(Sweep *sweep,float srange,float azim,
 278                                    float limit)
 279    {
 280    float  ccw_db_value,cw_db_value;
 281    double ccw_value,cw_value,value;
 282    double delta_angle;
 283    Ray    *ccw_ray,*cw_ray;
 284
 285    get_surrounding_ray(&ccw_ray,&cw_ray,sweep,azim);
 286
 287    /* Assume that ccw_ray and cw_ray will be non_NULL */
 288
 289    if((azim - ccw_ray->h.azimuth) > limit)
 290       {
 291       ccw_value = -1;
 292       }
 293    else
 294       {
 295       ccw_db_value = RSL_get_value_from_ray(ccw_ray,srange);
 296
 297       /* Check to make sure this is a valid value */
 298       if (ccw_db_value == BADVAL)
 299          {
 300          ccw_value = 0;
 301          }
 302       else
 303          {
 304          ccw_value = from_dB(ccw_db_value);
 305          }
 306       }
 307
 308    if((cw_ray->h.azimuth - azim) > limit)
 309       {
 310       cw_value = -1;
 311       }
 312    else
 313       {
 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)
 317          {
 318          cw_value = 0;
 319          }
 320       else
 321          {
 322          cw_value = from_dB(cw_db_value);
 323          }
 324       }
 325
 326    if((cw_value != -1) && (ccw_value != -1))
 327       {
 328       /* Both the clockwise ray and the counterclockwise
 329        * ray is valid.
 330        */
 331       delta_angle = angle_diff(ccw_ray->h.azimuth,cw_ray->h.azimuth);
 332
 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);
 335       }
 336    else if((cw_value == -1) && (ccw_value == -1))
 337       {
 338       /* Neither ray is valid. */
 339       value = -1;
 340       }
 341    else if(cw_value != -1)
 342       {
 343       /* cw_ray is only ray that is within limit. */
 344       value = cw_value;
 345       }
 346    else
 347       {
 348       /* ccw_ray is only ray that is within limit. */
 349       value = ccw_value;
 350       }
 351
 352    return value;
 353    }
 354
 355 /***************************************************/
 356 /*                                                 */
 357 /*          RSL_get_linear_value                   */
 358 /*                                                 */
 359 /*                                                 */
 360 /***************************************************/
 361 float RSL_get_linear_value(Volume *v,float srange,float azim,
 362                            float elev,float limit)
 363    {
 364    /* Compute bilinear value from four surrounding values
 365     * in Volume v.  The four surrounding values
 366     * are at a constant range.
 367     *
 368     * Limit is an angle used only to reject values
 369     * in the azimuth plane.
 370         */
 371    float  db_value;
 372    double value = 0;
 373    double up_value, down_value;
 374    double delta_angle;
 375    Sweep  *up_sweep,*down_sweep;
 376
 377    get_surrounding_sweep(&down_sweep,&up_sweep,v,elev);
 378
 379    /* Calculate interpolated value in sweep above
 380         * requested point.
 381         */
 382    if(up_sweep == NULL)
 383       {
 384       up_value = -1;
 385       }
 386    else
 387       {
 388       up_value = get_linear_value_from_sweep(up_sweep,srange,azim,limit);
 389       }
 390
 391    /* Calculate interpolated value in sweep below requested point.
 392         */
 393    if(down_sweep == NULL)
 394       {
 395       down_value = -1;
 396       }
 397    else
 398       {
 399       down_value = get_linear_value_from_sweep(down_sweep,srange,azim,limit);
 400       }
 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.
 404     */
 405    if((up_value != -1) && (down_value != -1))
 406       {
 407       delta_angle = angle_diff(up_sweep->h.elev,down_sweep->h.elev);
 408
 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);
 411       }
 412    else if((up_value == -1) && (down_value == -1))
 413       {
 414       value = -1;
 415       }
 416    else if(up_value != -1)
 417       {
 418       value = up_value;
 419       }
 420    else
 421       {
 422       value = down_value;
 423       }
 424
 425    /* Convert back to dB value and return. */
 426    if(value > 0)
 427       {
 428       db_value = (float)to_dB(value);
 429       return db_value;
 430       }
 431    else
 432       {
 433       return BADVAL;
 434       }
 435    }
 436