Pileus Git - ~andy/rsl/blob - src/rapic_routines.c

   1 #include <stdio.h>
   2 #include "rapic_routines.h"
   3 #include <string.h>
   4 #ifndef WIN32
   5 #include <netinet/in.h>
   6 #else
   7 #include "win32compat.h"
   8 #endif
   9
  10 void rapic_decode(unsigned char *inbuf, int inbytes, unsigned char *outbuf, int *outbytes,
  11                                   float *azim, float *elev, int *delta_time)
  12 {
  13   /* Decode RLE rapic buffer.
  14    *
  15    * Output to 'outbuf' and indicate the size w/ 'nout'.
  16    */
  17
  18   /* There is some risidual parsing to do:
  19    *
  20  *      @               <data>                  *
  21  *      <data>          AAA.A,EEE.E,TTT=LEN16D1D1D1NLD1D1D1NLNL
  22  *
  23  *      AAA.A           Azimuth in degrees
  24  *      EEE.E           Elevation in degress
  25  *      TTT             Delta time in seconds since the start of this scan
  26  *      LEN16           2 byte long length of radial
  27  *
  28  *      D1              Run length coded Data.
  29  *      NL              Null The Next Byte is count of NULL data.
  30  *      NLNL            End of this Radial
  31  *
  32  *         eg.          There will be no white space in the actual radial data.
  33  *
  34  *                      @066.1,010.6,004=002082B2817F8C84830048D72D0038
  35  *                                       999C0036202D35FD2C00238A99008AFE920000
  36  *                      Azimuth = 66.1
  37  *                      Elev    = 10.6
  38  *                      Dt      = 4 sec since the start
  39  *                      0020    = Bytes to follow
  40  *                      Data    = 82,B2,81,7F,8C,84,83
  41  *                      0048    = 48H null bytes
  42  *                      Data    = D7,2D
  43  *                      0038    = 38H null bytes
  44  *                      Data    = 99,9C
  45  *                      0036    = 36H Null bytes
  46  *                      ........................
  47  *                      0000    = End of Data.
  48  *      In  versions before 10.1
  49  *      @               <data>
  50  *      <data>          AAALEN16D1D1D1NLD1D1D1NLNL
  51  *
  52  *      AAA             Azimuth in degrees
  53  *      LEN16           2 byte long length of radial
  54  *
  55  *      D1              Run length coded Data.
  56  *      NL              Null The Next Byte is count of NULL data.
  57  *      NLNL            End of this Radial
  58  *
  59  *         eg.          There will be no white space in the actual radial data.
  60  *
  61  *                      @066002082B2817F8C84830048D72D0038
  62  *                          999C0036202D35FD2C00238A99008AFE920000
  63  *                      Azimuth = 66
  64  *                      0020    = Bytes to follow
  65  *                      Data    = 82,B2,81,7F,8C,84,83
  66  *                      0048    = 48H null bytes
  67  *                      Data    = D7,2D
  68  *                      0038    = 38H null bytes
  69  *                      Data    = 99,9C
  70  *                      0036    = 36H Null bytes
  71  *                      ........................
  72  *                      0000    = End of Data.
  73  *
  74  */
  75
  76
  77   /* The parser won't give us a line w/o '@' at the begining nor '\0\0'
  78    * at the end.  So we can be sure of that.
  79    */
  80
  81   int i;
  82   char prefix[16];
  83   unsigned short nnulls;
  84   short i16;
  85
  86   /* Find the '=' and start RLE decode from there. */
  87   *outbytes = 0;
  88   memset(prefix, '\0', sizeof(prefix));
  89   memcpy(prefix, &inbuf[1], 15);
  90
  91   sscanf(prefix, "%f,%f,%d", azim, elev, delta_time);
  92   /*  fprintf(stderr, "AZIM=%f, ELEV=%f, TTT=%d\n", *azim, *elev, *delta_time); */
  93
  94   /* Now, decode RLE. Don't care about 17,18 (they are the RLE buf size) */
  95   memcpy(&i16, &inbuf[17], 2);
  96   i16 = ntohs(i16);
  97   /*  fprintf(stderr, "Expecting %d bins\n", (int)i16); */
  98   i = 19;
  99   while (i<inbytes-2) {
 100         /*      fprintf(stderr, "i=%d byte=%2.2x(next %2.2x%2.2x) outbytes=%d\n", i, (unsigned char)inbuf[i], (unsigned char)inbuf[i+1], (unsigned char)inbuf[i+2], *outbytes); */
 101         if (inbuf[i] == '\0') { /* Next byte is a count of NULL's */
 102           i++;
 103           nnulls = (int)inbuf[i];
 104           /* fprintf(stderr, "NULL .. number of nulls=%4.4x\n", nnulls); */
 105           memset(&outbuf[*outbytes], '\0', (int)nnulls);
 106           *outbytes += nnulls;
 107         } else {
 108           /* fprintf(stderr, "Data\n"); */
 109           /* Data. */
 110           outbuf[*outbytes] = inbuf[i];
 111           *outbytes += 1;
 112         }
 113         i++;
 114   }
 115
 116   /* fprintf(stderr, "Decoded RLE: len=%d buf=<", *outbytes); */
 117   /* binprint(outbuf, *outbytes); */
 118   /* fprintf(stderr, ">\n"); */
 119 }
 120
 121
 122
 123 /*---------------------------------------------------------*/
 124 /*                                                         */
 125 /*                    binprint                             */
 126 /*                                                         */
 127 /*---------------------------------------------------------*/
 128 void binprint(char *s, int n)
 129 {
 130 int i;
 131
 132 for (i=0; i<n; i++)
 133          fprintf(stderr,"%c", s[i]);
 134
 135 }
 136
 137 /* Shut the damn linker up for librsl.so on Linux. This reference is
 138  * in the HDF library; don't need it, don't care.
 139  */
 140 void i_len(void)
 141 {
 142 }
 143
 144
 145 #include <time.h>
 146 void rapic_fix_time (Ray *ray)
 147 {
 148   struct tm the_time;
 149   float fsec;
 150   /* Fixes possible overflow values in month, day, year, hh, mm, ss */
 151   /* Normally, ss should be the overflow.  This code ensures end of
 152    * month, year and century are handled correctly by using the Unix
 153    * time functions.
 154    */
 155   if (ray == NULL) return;
 156   memset(&the_time, 0, sizeof(struct tm));
 157   the_time.tm_sec = ray->h.sec;
 158   fsec = ray->h.sec - the_time.tm_sec;
 159   the_time.tm_min  = ray->h.minute;
 160   the_time.tm_hour = ray->h.hour;
 161   the_time.tm_mon  = ray->h.month - 1;
 162   the_time.tm_year = ray->h.year - 1900;
 163   the_time.tm_mday = ray->h.day;
 164   the_time.tm_isdst = -1;
 165   (void) mktime(&the_time);
 166   /* The time is fixed. */
 167   ray->h.sec    = the_time.tm_sec;
 168   ray->h.sec   += fsec;
 169   ray->h.minute = the_time.tm_min;
 170   ray->h.hour   = the_time.tm_hour;
 171   ray->h.month  = the_time.tm_mon + 1;
 172   ray->h.year   = the_time.tm_year + 1900;
 173   ray->h.day    = the_time.tm_mday;
 174   return;
 175 }
 176
 177 void rapic_load_ray_header(Rapic_sweep_header rh, int iray, int isweep, float elev, float azim, Ray_header *h)
 178 {
 179   sscanf(rh.yyyymoddhhmmss,"%4d%2d%2d%2d%2d%2f",
 180                  &h->year,&h->month,&h->day,
 181                  &h->hour, &h->minute, &h->sec);
 182
 183   h->unam_rng = rh.end_range/1000.0;  /* Unambiguous range. (KM). */
 184   h->azimuth  = azim;   /* Azimuth angle. (degrees). Must be positive
 185                                         * 0=North, 90=east, -90/270=west.
 186                     * This angle is the mean azimuth for the whole ray.
 187                                         * Eg. for NSIG the beginning and end azimuths are
 188                                         *     averaged.
 189                                         */
 190   h->ray_num  = iray;     /* Ray no. within elevation scan. */
 191   h->elev     = rh.elev;  /* Elevation angle. (degrees). */
 192   h->elev_num = isweep;   /* Elevation no. within volume scan. */
 193
 194   h->range_bin1 = rh.start_range;       /* Range to first gate.(meters) */
 195   h->gate_size  = rh.range_resolution;  /* Data gate size (meters)*/
 196
 197   h->vel_res    = rh.range_resolution;   /* Doppler velocity resolution */
 198   h->sweep_rate = rh.anglerate/6.0;      /* Sweep rate. Full sweeps/min. */
 199
 200   h->prf       = rh.prf;          /* Pulse repetition frequency, in Hz. */
 201   h->azim_rate = rh.anglerate;    /* degrees/sec */
 202   h->fix_angle = elev;
 203   h->pitch     = 0;       /* Pitch angle. */
 204   h->roll      = 0;       /* Roll  angle. */
 205   h->heading   = 0;       /* Heading. */
 206   h->pitch_rate   = 0;            /* (angle/sec) */
 207   h->roll_rate    = 0;            /* (angle/sec) */
 208   h->heading_rate = 0;            /* (angle/sec) */
 209   h->lat          = rh.lat;       /* Latitude (degrees) */
 210   h->lon          = rh.lon;       /* Longitude (degrees) */
 211   h->alt          = rh.height;    /* Altitude (m) */
 212   h->rvc          = 0;            /* Radial velocity correction (m/sec) */
 213   h->vel_east     = 0;            /* Platform velocity to the east  (m/sec) */
 214   h->vel_north    = 0;            /* Platform velocity to the north (m/sec) */
 215   h->vel_up       = 0;            /* Platform velocity toward up    (m/sec) */
 216   h->pulse_count  = 0;            /* Pulses used in a single dwell time. */
 217   h->pulse_width  = rh.pulselen;  /* Pulse width (micro-sec). */
 218   h->beam_width   = rh.angle_resolution; /* Beamwidth in degrees. */
 219   h->frequency    = rh.freq/1000.0; /* Carrier freq. GHz. */
 220   h->wavelength   = 0;            /* Wavelength. Meters. */
 221   h->nyq_vel      = rh.nyquist;   /* Nyquist velocity. m/s */
 222
 223   return;
 224 }
 225
 226 extern float rapic_nyquist;
 227
 228 float RAPIC_DZ_F(unsigned char x) {
 229   if (x == 0) return NOECHO;
 230   return (((float)x-64)/2.0);  /* rapic -> float */
 231 }
 232
 233 float RAPIC_VR_F(unsigned char x) {
 234   if (x == 0) return BADVAL;
 235   return (((float)((int)x-128))/128.0*rapic_nyquist);  /* rapic -> float */
 236 }
 237
 238 float RAPIC_SW_F(unsigned char x) {
 239   if (x == 0) return NOECHO;
 240   return (((float)x)/256.0*rapic_nyquist);  /* rapic -> float */
 241 }
 242
 243 float RAPIC_ZT_F(unsigned char x) {
 244   return RAPIC_DZ_F(x);
 245 }
 246
 247 float RAPIC_ZD_F(unsigned char x) {
 248   if (x == 0) return NOECHO;
 249   return (((float)x-128)/16.0);  /* rapic -> float */
 250 }
 251
 252 /* USE RSL INDEXING! */
 253 static  float (*RAPIC_f_list[])(unsigned char x) = {RAPIC_DZ_F,
 254                                                                                                         RAPIC_VR_F,
 255                                                                                                         RAPIC_SW_F,
 256                                                                                                         NULL,
 257                                                                                                         RAPIC_ZT_F,
 258                                                                                                         NULL,
 259                                                                                                         NULL,
 260                                                                                                         RAPIC_ZD_F};
 261
 262
 263 void rapic_load_ray_data(unsigned char *buf, int bufsize, int ifield, Ray *ray)
 264 {
 265   /* ifield is the RSL numbering scheme for field types.  The conversion
 266    * is done in rapic.y.  In other words, we've already converted rapic
 267    * index to RSL indexes.
 268    */
 269   int i;
 270   for (i=0; i<bufsize; i++) {
 271         ray->range[i] = ray->h.invf(RAPIC_f_list[ifield](buf[i]));
 272         /*      fprintf(stderr,"i=%d ifield=%d, buf[%d]=%2.2x, ray->range[%d]=%4.4x value=%f\n", i,ifield,i,buf[i],i,ray->range[i], RAPIC_f_list[ifield](buf[i]) ); */
 273   }
 274   ray->h.nbins = bufsize;
 275 }
 276
 277 Radar *fill_header(Radar *radar)
 278 {
 279   /* Learn as much as possible from the Ray headers.  Place this
 280    * information into radar->h.xxxxxx
 281    */
 282   Ray *ray;
 283   Volume *volume;
 284   int i;
 285   float tmp;
 286
 287   volume = NULL;
 288   if (radar == NULL) return NULL;
 289   for (i=0; i<radar->h.nvolumes && !(volume = radar->v[i]); i++)
 290         ;
 291   if (volume == NULL) return NULL;
 292
 293   ray = RSL_get_first_ray_of_volume(volume);
 294   if (ray  == NULL) return NULL;
 295
 296   radar->h.month = ray->h.month;
 297   radar->h.day   = ray->h.day;
 298   radar->h.year  = ray->h.year;
 299   radar->h.hour  = ray->h.hour;
 300   radar->h.minute= ray->h.minute;
 301   radar->h.sec   = ray->h.sec; /* Second plus fractional part. */
 302   sprintf(radar->h.radar_type, "rapic"); /* Type of radar. */
 303                   /* nvolumes   is already filled in YACC. */
 304                   /* number     is already filled in YACC. */
 305                   /* name       is already filled in YACC. */
 306                   /* radar_name is already filled in YACC. */
 307   /*  radar->h.city[15];     */ /* Not available from RAPIC. */
 308   /*  radar->h.state[3];     */ /* Not available from RAPIC. */
 309   /*  radar->h.country[15];  */ /* Not available from RAPIC. */
 310
 311    /** Latitude deg, min, sec **/
 312   radar->h.latd = (int)ray->h.lat;
 313   tmp = (ray->h.lat - radar->h.latd) * 60.0;
 314   radar->h.latm = (int)tmp;
 315   radar->h.lats = (int)((tmp - radar->h.latm) * 60.0);
 316    /** Longitude deg, min, sec **/
 317   radar->h.lond = (int)ray->h.lon;
 318   tmp = (ray->h.lon - radar->h.lond) * 60.0;
 319   radar->h.lonm = (int)tmp;
 320   radar->h.lons = (int)((tmp - radar->h.lonm) * 60.0);
 321
 322   radar->h.height = ray->h.alt; /* height of site in meters above sea level*/
 323   radar->h.spulse = 0; /* length of short pulse (ns)*/
 324   radar->h.lpulse = 0; /* length of long pulse (ns) */
 325
 326   return radar;
 327 }