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

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             John H. Merritt
   6             Space Applications Corporation
   7             Vienna, Virginia
   8
   9     This library is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU Library General Public
  11     License as published by the Free Software Foundation; either
  12     version 2 of the License, or (at your option) any later version.
  13
  14     This library is distributed in the hope that it will be useful,
  15     but WITHOUT ANY WARRANTY; without even the implied warranty of
  16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17     Library General Public License for more details.
  18
  19     You should have received a copy of the GNU Library General Public
  20     License along with this library; if not, write to the Free
  21     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  22 */
  23 #include <stdlib.h>
  24 #include <math.h>
  25
  26 #  define M_PI          3.14159265358979323846
  27
  28 #include "rsl.h"
  29 /**********************************************************************/
  30 /*                                                                    */
  31 /*                 RSL_area_of_ray                                    */
  32 /*                 RSL_fractional_area_of_sweep                       */
  33 /*                                                                    */
  34 /*    Compute fractional area of a sweep given a range of dBZ and     */
  35 /*    a maximum range in km.                                          */
  36 /*                                                                    */
  37 /*  By: John Merritt                                                  */
  38 /*      Space Applications Corporation                                */
  39 /*      April 7, 1994                                                 */
  40 /**********************************************************************/
  41
  42 float get_pixel_area(Ray *r, float min_range, float max_range)
  43 {
  44   float h1, h2, r1, r2;  /* height and radius in km*/
  45   float volume;
  46   float theta;  /* in radian */
  47
  48   /* returns volume area of ray between min range and max range
  49    */
  50   if (r == NULL) return 0.0;
  51
  52   h1 = min_range;
  53   h2 = max_range;
  54   /* convert from degree to radian */
  55   theta = (r->h.beam_width/2) * 2 * M_PI / 360;
  56
  57   r1 = h1 * tan(theta);
  58   r2 = h2 * tan(theta);
  59   volume = (M_PI/3) * (pow((double)r2, (double) 2) * h2 - pow((double) r1, (double) 2) *h1);
  60
  61   return volume;
  62 } /* get_pixel_area */
  63
  64
  65 float RSL_area_of_ray(Ray *r, float lo, float hi, float min_range, float max_range)
  66 {
  67   int i;
  68   float start_km;
  69   float xdBZ;
  70   float r1, r2, area;
  71   int nbins, bin1;
  72   float binsize;
  73
  74   if (r == NULL) return 0.0;
  75   /* Check if min_range is closer to the radar than the first bin.
  76    * If so, we can't use it.
  77    */
  78
  79   if (hi <= 0.0)
  80         hi = 100.0;  /* default to max dbz*/
  81   start_km = r->h.range_bin1/1000.0;
  82   binsize = r->h.gate_size/1000.0;
  83
  84   nbins = ((max_range - start_km) / binsize);
  85   bin1 = (min_range - start_km) / binsize;
  86
  87   if (bin1 < 0) bin1 = 0;
  88   if (nbins > r->h.nbins) nbins = r->h.nbins;
  89
  90   /* Compute the number of pixels with lo < dBZ <= hi */
  91   area = 0.0;
  92
  93   for(i=bin1; i<nbins-1; i++) {
  94         xdBZ = r->h.f(r->range[i]);
  95         if(lo < xdBZ && xdBZ <= hi) { /*MAX_DBZ = hi (typically 70?) */
  96           /* get r1 and r2 in km */
  97           r1 = i * binsize + start_km;
  98           r2 = (i+1) * binsize + start_km;
  99           area += get_pixel_area(r, r1, r2);
 100         }
 101         r1 = i * binsize + start_km;
 102   }
 103
 104   return area;
 105 }
 106
 107
 108
 109 float RSL_fractional_area_of_sweep(Sweep *s, float lo, float hi, float min_rng, float max_rng)
 110 {
 111 /*
 112  *  Compute the fractional area of the Sweep.
 113  *
 114  *  This doesn't take care of a sector (window).  The caller will
 115  *  have to multiply the answer by some appropriate constant.
 116  *
 117  */
 118
 119
 120   float sweep_area;
 121   float area;
 122   float frac_area;
 123   float total_sys_area=0.0;
 124   int iazm;
 125
 126   if (s == NULL) return 0.0;
 127
 128   sweep_area = 0;
 129   for (iazm=0; iazm<s->h.nrays; iazm++) {
 130         /* get total system volume area */
 131         total_sys_area += RSL_area_of_ray(s->ray[iazm], -71, 0.0, min_rng, max_rng);
 132         area = RSL_area_of_ray(s->ray[iazm], lo, hi, min_rng, max_rng);
 133         sweep_area += area;
 134 /*
 135         fprintf(stderr,"iazm = %d, totalsysarea = %f, area = %f, sweep_area = %f\n",
 136                    iazm, total_sys_area, area, sweep_area);
 137 */
 138   }
 139   frac_area = sweep_area / total_sys_area;
 140   return (float) frac_area;
 141 }
 142