Ticket 742: Remove old prototype code. hudson-beamb-28-SUCCESS
authorAnders Henja <anders@baltrad.eu>
Wed, 4 Jan 2012 15:31:15 +0000 (16:31 +0100)
committerAnders Henja <anders@baltrad.eu>
Wed, 4 Jan 2012 15:31:15 +0000 (16:31 +0100)
lib/beamb.c [deleted file]
lib/beamb_fn.c [deleted file]
lib/beamb_fn.h [deleted file]
lib/beamb_map.c [deleted file]
lib/beamb_map.h [deleted file]

diff --git a/lib/beamb.c b/lib/beamb.c
deleted file mode 100644 (file)
index b8af09a..0000000
+++ /dev/null
@@ -1,344 +0,0 @@
-/* --------------------------------------------------------------------
-Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
-
-This is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-\13
-This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License along with HLHDF.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------*/
-/**
- * Beam-blockage analysis
- * @file
- * @author Lars Norin (Swedish Meteorological and Hydrological Institute, SMHI)
- * @date 2011-10-01
- */
-#include <math.h>
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-
-#ifndef M_PI
-/**
- * Value of PI. Can be replaced with that found in PROJ.4 ...
- */
-#define M_PI 3.14159265358979323846
-#endif
-
-/**
- * Struct containing map parameters.
- * ulxmap - longitude of the center of the upper-left pixel (decimal degrees)
- * ulymap - latitude  of the center of the upper-left pixel (decimal degrees)
- * nbits - number of bits per pixel (16 for a DEM)
- * nrows - number of rows in the image
- * ncols - number of columns in the image
- * xdim - x dimension of a pixel in geographic units (decimal degrees)
- * ydim - y dimension of a pixel in geographic units (decimal degrees)
- */
-typedef struct
-{
-    double ulxmap;
-    double ulymap;
-    int nbits;
-    int nrows;
-    int ncols;
-    double xdim;
-    double ydim;
-} map_info;
-
-#include "beamb_fn.h"
-#include "beamb_map.h"
-
-/**
- * Calculates percentage of beam power blocked.
- * Result is written to file.
- * @param [in] hdf5file - hdf5 file containing metadata on radar
- * @param [in] dBlim - limit of Gaussian approximation of main lobe
- */
-int main(void)
-{
-    /* Declaration of variables */
-    char *radar_name, filedir[100], filedir1[100], filedir2[100];
-    int filenr;
-    double radar[3];
-    double el, beamwidth, phi0, az_step;
-    double latRadar, lonRadar, hRadar;
-    double r_step, r_min, r_max, dBlim;
-    double range_len, a_len;
-    int ri, ai;
-
-    /* Directory where map files are located */
-    strcpy(filedir,"/data/proj/radar/lnorin/baltrad/beam_blockage/gtopo30/");
-    strcpy(filedir1,"/data/proj/radar/lnorin/baltrad/beam_blockage/gtopo30/");
-    strcpy(filedir2,"/data/proj/radar/lnorin/baltrad/beam_blockage/gtopo30/");
-
-    /* Limit of Gaussian approximation of main lobe (given as parameter to main?) */
-    dBlim=-20.;
-    
-    /* Name of the radar (should be read from hdf5 file) */
-/*    radar_name = "arl";*/
-    radar_name = "var";
-    
-    /* Properties of the radar (should be read from hdf5 file) */
-    if (strcmp(radar_name,"arl") == 0)
-    {
-        radar[0] = 59.654437083; /* Latitude [degrees] */
-        radar[1] = 17.946310106; /* Longitude [degrees] */
-        radar[2] = 73.51;        /* Altitude [m] */
-    }
-    else if (strcmp(radar_name,"var") == 0)
-    {
-        radar[0] = 58.255645047; /* Latitude [degrees] */
-        radar[1] = 12.826024108; /* Longitude [degrees] */
-        radar[2] = 163.61;       /* Altitude [m] */
-    }
-    else if (strcmp(radar_name,"vil") == 0)
-    {
-      radar[0] = 58.1059;      /* Latitude [degrees] */
-      radar[1] = 15.9363;      /* Longitude [degrees] */
-      radar[2] = 223.0;        /* Altitude [m] */
-    }
-    else
-    {
-        printf("Unknown radar. Exiting.");
-        return 0;
-    }
-
-    /* Properties of the radar (should be read from hdf5 file) */
-    if (strcmp(radar_name,"hur") == 0)
-    {
-            /*** Norwegian radars ***/
-            el = .5;          /* Elevation angles [degrees] */
-            beamwidth = 1.;   /* Beamwidth [degrees] */
-            phi0 = 90.;       /* Start position for sweep [degrees] */
-            az_step = 1.;
-    }
-    else
-    {
-            /*** Swedish radars ***/
-            el = .5;          /* Elevation angles [degrees] */
-            beamwidth = .9;   /* Beamwidth [degrees] */
-            phi0 = 90.;       /* Start position for sweep [degrees] */
-            az_step = 360./420.;
-    }
-
-    /* Coordinates of the radar */
-    latRadar = radar[0];
-    lonRadar = radar[1];
-    hRadar = radar[2];
-    
-    /* Range resolution [m] (should be read from hdf5 file) */
-    r_step = 2.e3;
-    r_min = 0.;
-    r_max = 240.e3;
-
-    /* Number of range bins and azimuth gates */
-    range_len = floor((r_max-r_min)/r_step);
-    ri = (int) range_len;
-    a_len = 360./az_step;
-    ai = (int) a_len;
-    
-    /* Allocating variables */
-    double *range = malloc(sizeof(double)*ri);
-    double *azimuth = malloc(sizeof(double)*ai);
-    
-    /* Polar grid of the radar, centered in each pixel */
-    int i;
-    for (i=0; i<ai; i++)
-        azimuth[i] = az_step/2 + i*az_step;
-    
-    for (i=0; i<ri; i++)
-        range[i] = r_min + r_step/2 + i*r_step;
-    
-    /* Find out which map to read */
-    filenr = getTopo(latRadar, lonRadar, range[ri-1]);
-    
-    /* Allocating variables */
-    map_info *map = malloc(sizeof(*map));
-    map_info *map1 = malloc(sizeof(*map1));
-    map_info *map2 = malloc(sizeof(*map2));
-    
-    int nrows, ncols, nbytes;
-    double ulxmap, ulymap, xdim, ydim;
-    
-    switch (filenr)
-    {
-        case 1:
-            /* One map covers radar grid */
-            strcat(filedir,"W020N90"); 
-            /* Read map header */
-            readHdr(map, filedir);
-            nrows = map->nrows;
-            ncols = map->ncols;
-            nbytes = (map->nbits)/8;
-            xdim = map->xdim;
-            ydim = map->ydim;
-            ulxmap = map->ulxmap; 
-            ulymap = map->ulymap;
-            break;
-        case 2:
-            /* One map covers radar grid */
-            strcat(filedir,"E020N90"); 
-            /* Read map header */
-            readHdr(map, filedir);
-            nrows = map->nrows;
-            ncols = map->ncols;
-            nbytes = (map->nbits)/8;
-            xdim = map->xdim;
-            ydim = map->ydim;
-            ulxmap = map->ulxmap; 
-            ulymap = map->ulymap;
-            break;
-        case 3:
-            /* Two maps must be used to cover radar grid */
-            strcat(filedir1,"W020N90"); 
-            strcat(filedir2,"E020N90"); 
-            /* Read map headers */
-            readHdr(map1, filedir1);
-            readHdr(map2, filedir2);
-            /* Check that number of rows and spacing are the same */
-            if(map1->nrows == map2->nrows && map1->xdim == map2-> xdim \
-                && map1->ydim == map2-> ydim && map1->nbits == map2->nbits)
-            {
-                nrows = map1->nrows;
-                ncols = map1->ncols + map2->ncols;
-                nbytes = (map1->nbits)/8;
-                xdim = map1->xdim;
-                ydim = map1->ydim;
-                ulxmap = map1->ulxmap; 
-                ulymap = map1->ulymap;
-            }
-            break;
-    }
-    
-    /* Should maybe be extended to cover other format than int16 */
-    if(nbytes != 2)
-        printf("Warning! Unexpected file format.");
-    
-    /* Allocating variables */
-    short int *data = malloc(nrows*ncols*nbytes);
-    
-    /* Read map(s) */
-    switch (filenr)
-    {
-        case 1:
-            readMap(data, map, filedir);
-            break;
-        case 2:
-            readMap(data, map, filedir);
-            break;
-        case 3:
-            readMap2(data, map1, map2, filedir1, filedir2);
-            break;
-    }
-    
-    /* Allocating variables */
-    double *ground_range = malloc(sizeof(double)*ri);
-    double *lat = malloc(sizeof(double)*ri*ai);
-    double *lon = malloc(sizeof(double)*ri*ai);
-    double *zi = malloc(sizeof(double)*ri*ai);
-    
-    /* Find the projected range on the ground [m] */
-    compute_ground_range(ground_range, range, latRadar, hRadar, el, ri);
-    
-    /* Convert radar polar grid to lat/long */
-    polar2latlon(lat, lon, latRadar, lonRadar, ground_range, azimuth, ri, ai);
-
-    /* Find indices of lat and lon limits of radar grid */
-    int lon_min, lon_max, lat_min, lat_max, lon1, lon2, lat1, lat2;
-    
-    lon_min = floor( (min_double(lon,ri*ai)-ulxmap) / xdim ) - 1;
-    lon_max = ceil( (max_double(lon,ri*ai)-ulxmap) / xdim ) + 1;
-    lat_max = floor( (ulymap-max_double(lat,ri*ai)) / ydim ) - 1;
-    lat_min = ceil( (ulymap-min_double(lat,ri*ai)) / ydim ) + 1;
-
-    /* Check that we are not outside map */
-    lon1 = lon_min > 1 ? lon_min : 1;
-    lon2 = lon_max < ncols ? lon_max : ncols;
-    lat1 = lat_max > 1 ? lat_max : 1;
-    lat2 = lat_min < nrows ? lat_min : nrows;
-    
-    /* Declaration of variables */
-    int nlon = lon2 - lon1 +1;
-    int nlat = lat2 - lat1 +1;
-
-    /* Allocating variables */
-    double *data_small = malloc(sizeof(double)*nlon*nlat);
-    
-    /* Cut out smaller part of map */
-    int j, k = 0;
-    for(i = lat1-1; i < lat2; i++)
-    {
-        for(j = lon1-1; j < lon2; j++)
-        {
-            if(data[i*ncols+j] < 0)
-            {
-                data_small[k] = 0;
-                k++;
-            }
-            else
-            {
-                data_small[k] = (double) data[i*ncols+j];
-                k++;
-            }
-        }
-    }
-    
-   /* Allocating variables */
-    double *lon_map = malloc(sizeof(double)*nlon);
-    double *lat_map = malloc(sizeof(double)*nlat);
-    
-   /* Define the grid for the small map */
-    for(i=0; i<nlon; i++)
-    {
-        lon_map[i] = ulxmap+(lon1-1)*xdim + i*xdim;
-    }
-    for(i=0; i<nlat; i++)
-    {
-        lat_map[i] = ulymap-(lat1-1)*ydim - i*ydim;
-    }
-
-    /* Interpolate map from original grid to radar grid. */
-    BilinearInterpolation(lon_map, lat_map, data_small, nlon, nlat, \
-                         lon, lat, zi, ri, ai);
-   
-    /* Allocating variables */
-    double *bb = malloc(sizeof(double)*ri*ai);
-    
-    /* Compute blockage */
-    getBlockage(bb, zi, latRadar, hRadar, ground_range, el, beamwidth, \
-                dBlim, ri, ai);
-    
-    /* Write to file */
-    FILE *fp1;
-    int n;
-    if((fp1 = fopen("/data/proj/radar/lnorin/baltrad/beam_blockage/bb.dat", "w")))
-    {
-        n = fwrite(bb, sizeof(double), ri*ai, fp1);
-        fclose(fp1);
-    }
-    else
-    {
-        printf("Error writing file.");
-    }
-    
-    /* Free memory */
-    free(range);
-    free(ground_range);
-    free(azimuth);
-    free(lat);
-    free(lon);
-    free(map);
-    free(map1);
-    free(map2);
-    free(data);
-    free(data_small);
-    free(lon_map);
-    free(lat_map);
-    free(zi);
-    free(bb);
-    
-    return 0;
-}
-
-
diff --git a/lib/beamb_fn.c b/lib/beamb_fn.c
deleted file mode 100644 (file)
index 2cc6d64..0000000
+++ /dev/null
@@ -1,247 +0,0 @@
-/* --------------------------------------------------------------------
-Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
-
-This is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-
-This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License along with HLHDF.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------*/
-/**
- * Convenience functions for beam-blockage analysis
- * @file
- * @author Lars Norin (Swedish Meteorological and Hydrological Institute, SMHI)
- * @date 2011-10-01
- */
-#include <math.h>
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-
-#ifndef M_PI
-/**
- * Value of PI. Can be replaced with that found in PROJ.4 ...
- */
-#define M_PI 3.14159265358979323846
-#endif
-
-/**
- * Converts degrees to radians
- * @param[in] deg - input value expressed in degrees
- * @returns value converted to radians
- */
-double deg2rad(double deg)
-{
-    return deg*M_PI/180.;
-}
-
-/**
- * Convert radians to degrees
- * @param[in] rad - input value expressed in radians
- * @returns value converted to degrees
- */
-double rad2deg(double rad)
-{
-    return rad*180./M_PI;
-}
-
-/**
- * Get simplified Earth radius for a certain latitude
- * @param[in] lat0 - input value expressed in degrees
- * @returns Earth radius at lat0 in meters
- */
-double getEarthRadius(double lat0)
-{
-    double R_EQU, R_POL, a, b, radius;
-    R_EQU = 6378160.; /* [m] */
-    R_POL = 6356780.; /* [m] */
-    a = sin(deg2rad(lat0))*R_POL;
-    b = cos(deg2rad(lat0))*R_EQU;
-    radius = sqrt(pow(a,2)+pow(b,2)); /* [m] */
-    return radius;
-}
-
-/**
- * Get range from radar, projected on surface
- * @param[in] *ground_range - projected range in meters
- * @param[in] *range - radar range in meters
- * @param[in] lat0 - latitude in degrees
- * @param[in] alt0 - radar altitude in meters
- * @param[in] el - elevation of radar beam in degrees
- * @param[in] r_len - length of vector range
- */
-void compute_ground_range(double *ground_range, double *range, double lat0, \
-                         double alt0, double el, int r_len)
-{
-    double dndh = -3.9e-8, RE, R, A, H, gamma;
-    int i;
-    RE = getEarthRadius(lat0);
-    R = 1./((1./RE) + dndh);
-    A = R + alt0;
-    for(i = 0; i < r_len; i++)
-    {
-        H = sqrt(pow(A,2) + pow(range[i],2) + 2*A*range[i]*sin(deg2rad(el))) \
-            -A;
-        gamma = asin(range[i]*sin(M_PI/2+deg2rad(el))/(A+H));
-        ground_range[i] = A * gamma;
-    }
-    return;
-}
-
-/**
- * Given a start point, initial bearing, and distance, this will calculate the
- * destination point and final bearing travelling along a (shortest distance) great
- * circle arc. See http://www.movable-type.co.uk/scripts/latlong.html
- * @param[in] *lat - vector of latitudes in degrees
- * @param[in] *lon - vector of longitudes in degrees
- * @param[in] latR - radar latitude in degrees
- * @param[in] lonR - radar longitude in meters
- * @param[in] *range - elevation of radar beam in degrees
- * @param[in] *azimuth - length of vector range
- * @param[in] r_len - length of vector range
- * @param[in] a_len - length of vector range
- */
-void polar2latlon(double *lat, double *lon, double latR, double lonR, \
-                  double *range, double *azimuth, int r_len, int a_len)
-{
-    double RE;
-    int i, j;
-    RE=getEarthRadius(latR);
-    
-    for(i=0; i<a_len; i++)
-    {
-        for(j=0; j<r_len; j++)
-        {
-            lat[i*r_len+j] = rad2deg( asin( sin(deg2rad(latR)) * cos(range[j]/RE) + \
-                            cos(deg2rad(latR)) * sin(range[j]/RE) * \
-                            cos(deg2rad(azimuth[i])) ) );
-            lon[i*r_len+j] = lonR + rad2deg( atan2( cos(deg2rad(latR)) * \
-                            sin(range[j]/RE) * sin(deg2rad(azimuth[i])) , \
-                            cos(range[j]/RE) - sin(deg2rad(latR)) * \
-                            sin(deg2rad(lat[i*r_len+j])) ) );
-        }
-    }
-}
-
-/**
- * Find minimum value in vector
- * @param[in] *a - pointer to vector of values
- * @param[in] n - length of vector a
- * @returns minimum value of a
- */
-double min_double(double *a, int n)
-{
-    double b = a[0];
-    int i;
-    for(i=0; i<n; i++)
-    {
-        if (a[i] < b)
-        {
-            b = a[i];
-        }
-    }
-    return b;
-}
-
-/**
- * Find maximum value in vector
- * @param[in] *a - pointer to vector of values
- * @param[in] n - length of vector a
- * @returns maximum value of a
- */
-double max_double(double *a, int n)
-{
-    double b = a[0];
-    int i;
-    for(i=0; i<n; i++)
-    {
-        if (a[i] > b)
-        {
-            b = a[i];
-        }
-    }
-    return b;
-}
-
-/**
- * Bilinear interpolation, see e.g. 
- * http://en.wikipedia.org/wiki/Bilinear_interpolation
- * Given x1, y1, z1 and x2, y2 function finds z2
- * @param[in] *x1 - pointer to vector of abcissae
- * @param[in] *y1 - pointer to vector of ordinates
- * @param[in] *z1 - pointer to vector of values
- * @param[in] x1_max - length of vector x1
- * @param[in] y1_max - length of vector y1
- * @param[in] *x2 - pointer to vector of abcissae
- * @param[in] *y2 - pointer to vector of ordinates
- * @param[in] *z2 - pointer to vector of values
- * @param[in] x2_max - length of vector x2
- * @param[in] y2_max - length of vector y2
- */
-void BilinearInterpolation(double *x1, double *y1, double *z1, int x1_max, \
-                           int y1_max, double *x2, double *y2, double *z2, \
-                           int x2_max, int y2_max)
-{
-    int i, j, m1, m2, n1, n2;
-    double q1, q2, q3, q4, dx, dy, dxdy;
-    double x1min = min_double(x1,x1_max), x1max = max_double(x1,x1_max), \
-           y1min = min_double(y1,y1_max), y1max = max_double(y1,y1_max);
-    dx = x1[1]-x1[0];
-    dy = y1[1]-y1[0];
-    dxdy = dx*dy;
-    
-    for (j=0; j<y2_max; j++)
-    {
-        for (i=0; i<x2_max; i++)
-        {
-            /* If (x2,y2) is outside ([x_min:x_max],[y_min:y_max]),
-               let value be zero (no extrapolation) */
-            if (x2[j*x2_max+i]<x1min || x2[j*x2_max+i]>x1max || \
-                y2[j*x2_max+i]<y1min || y2[j*x2_max+i]>y1max)
-            {
-                z2[j*x2_max+i] = 0.0;
-                printf("Warning! Function BilinearInterpolation does not \
-                        perform extrapolation.");
-                continue;
-            }
-            
-            /* Find coordinates of surrounding grid points */
-            m1 = (int) floor( (x2[j*x2_max+i] - x1[0]) / dx );
-            n1 = (int) floor( (y2[j*x2_max+i] - y1[0]) / dy );
-            m2 = m1 + 1;
-            n2 = n1 + 1;
-            
-            q1 = (x1[m2]-x2[j*x2_max+i]) * (y1[n2]-y2[j*x2_max+i]) / dxdy;
-            q2 = (x2[j*x2_max+i]-x1[m1]) * (y1[n2]-y2[j*x2_max+i]) / dxdy;
-            q3 = (x1[m2]-x2[j*x2_max+i]) * (y2[j*x2_max+i]-y1[n1]) / dxdy;
-            q4 = (x2[j*x2_max+i]-x1[m1]) * (y2[j*x2_max+i]-y1[n1]) / dxdy;
-            
-            z2[j*x2_max+i] = z1[n1*x1_max+m1]*q1 + \
-                             z1[n2*x1_max+m1]*q3 + \
-                             z1[n1*x1_max+m2]*q2 + \
-                             z1[n2*x1_max+m2]*q4;
-        }
-    }
-    return;
-}
-
-/**
- * Running cumulative maximum
- * @param[in] *a - pointer to vector of values
- * @param[in] a_max - length of vector a
- */
-void cummax(double *a, int a_max)
-{
-    double t;
-    int i;
-    t = a[0];
-    for(i=1; i<a_max; i++)
-    {
-        if(a[i]<t)
-            a[i] = t;
-        else
-            t = a[i];
-    }
-    return;
-}
-
diff --git a/lib/beamb_fn.h b/lib/beamb_fn.h
deleted file mode 100644 (file)
index 9f97dc9..0000000
+++ /dev/null
@@ -1,24 +0,0 @@
-/* Various functions used by beamb.c */
-
-double deg2rad(double deg);
-
-double rad2deg(double rad);
-
-double getEarthRadius(double lat);
-
-void compute_ground_range(double *slant_range, double *range, double lat0, \
-                         double alt0, double el, int range_len);
-
-void polar2latlon(double *lat, double *lon, double latR, double lonR, \
-                  double *range, double *azimuth, int r_len, int a_len);
-
-double min_double(double *a, int n);
-
-double max_double(double *a, int n);
-
-void BilinearInterpolation(double *x1, double *y1, double *z1, int x1_max, \
-                           int y1_max, double *x2, double *y2, double *z2, \
-                           int x2_max, int y2_max);
-
-void cummax(double *a, int a_max);
-
diff --git a/lib/beamb_map.c b/lib/beamb_map.c
deleted file mode 100644 (file)
index ab2b472..0000000
+++ /dev/null
@@ -1,334 +0,0 @@
-/* --------------------------------------------------------------------
-Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
-
-This is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-
-This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License along with HLHDF.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------*/
-/**
- * Map reading-functions for beam-blockage analysis
- * @file
- * @author Lars Norin (Swedish Meteorological and Hydrological Institute, SMHI)
- * @date 2011-10-01
- */
-#include <math.h>
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <arpa/inet.h>
-
-/**
- * Struct containing map parameters.
- * ulxmap - longitude of the center of the upper-left pixel (decimal degrees)
- * ulymap - latitude  of the center of the upper-left pixel (decimal degrees)
- * nbits - number of bits per pixel (16 for a DEM)
- * nrows - number of rows in the image
- * ncols - number of columns in the image
- * xdim - x dimension of a pixel in geographic units (decimal degrees)
- * ydim - y dimension of a pixel in geographic units (decimal degrees)
- */
-typedef struct
-{
-    double ulxmap;
-    double ulymap;
-    int nbits;
-    int nrows;
-    int ncols;
-    double xdim;
-    double ydim;
-} map_info;
-
-#include "beamb_fn.h"
-
-#ifndef M_PI
-/**
- * Value of PI. Can be replaced with that found in PROJ.4 ...
- */
-#define M_PI 3.14159265358979323846
-#endif
-
-/**
- * Find out which maps are needed to cover given area
- * @param[in] lat - latitude of radar in degrees
- * @param[in] lon - longitude of radar in degrees
- * @param[in] d - maximum range of radar in meters
- * @returns flag corresponding to map to be read
- */
-int getTopo(double lat, double lon, double d)
-{
-    double latRad, lonRad, RE, lat_e, lon_e, lat_w, lon_w, lat_n, lat_s;
-    
-    lonRad = deg2rad(lon);
-    latRad = deg2rad(lat);
-    RE = getEarthRadius(lat);
-    
-    /* Find out which map is needed
-       see http://www.movable-type.co.uk/scripts/latlong.html */
-    lat_e = asin(sin(latRad) * cos(d/RE) + cos(latRad) * sin(d/RE) * cos(M_PI/2.));
-    lon_e = lonRad + atan2( sin(M_PI/2.) * sin(d/RE) * cos(latRad), \
-            cos(d/RE) - sin(latRad) * sin(lat_e) );
-    lat_w = asin( sin(latRad) * cos(d/RE) + cos(latRad) * sin(d/RE) * cos(3.*M_PI/2.) );
-    lon_w = lonRad + atan2( sin(3.*M_PI/2.) * sin(d/RE) * cos(latRad), \
-            cos(d/RE) - sin(latRad) * sin(lat_w) );
-    lat_n = asin( sin(latRad) * cos(d/RE) + cos(latRad) * sin(d/RE) * cos(0.) );
-    lat_s = asin( sin(latRad) * cos(d/RE) + cos(latRad) * sin(d/RE) * cos(M_PI) );
-    
-    /* Check if maps cover the area */
-    if(rad2deg(lat_n)>80. || rad2deg(lat_s)<40. || rad2deg(lon_w)<-20. || \
-       rad2deg(lon_e)>60.)
-    {
-        printf("Maps do not cover area. Exiting.\n");
-        return -1;
-    }
-    
-    if(rad2deg(lon_e)<=20.)
-    {
-        /* Read W020N90 */
-        return 1;
-    }
-    else if(rad2deg(lon_w)>20.)
-    {
-        /* Read E020N90 */
-        return 2;
-    }
-    else
-    {
-        /* Read both maps */
-        return 3;
-    }
-    printf("...done!\n");
-}
-
-/**
- * Read header of gtopo30 files
- * @param[in] *map - pointer to struct containing map information
- * @param[in] *filename - name of file to be read
- */
-void readHdr(map_info *map, char *filename)
-{
-    /* Declaration of variables */
-    char filename_copy[100];
-    char *s, line[100];
-    int n = 15, f_int;
-    double f;
-    
-    strcpy(filename_copy,filename);
-    FILE *fp = fopen(strcat(filename_copy,".HDR"),"rt");
-    if(fp != (FILE *)NULL)
-    {
-        /* Get number of rows and columns, coordinates, and lat/lon increments */
-        s = (char *)malloc(n);
-        while(fgets(line, sizeof(line), fp) > 0)
-        {
-            if(sscanf(line, "%s %lf", s, &f) == 2)
-            {
-                if(strcmp(s,"NROWS") == 0)
-                {
-                    f_int = (int) f;
-                    map->nrows = f_int;
-                }
-                else if(strcmp(s,"NCOLS") == 0)
-                {
-                    f_int = (int) f;
-                    map->ncols = f_int;
-                }
-                else if(strcmp(s,"NBITS") == 0)
-                {
-                    f_int = (int) f;
-                    map->nbits = f_int;
-                }
-                else if(strcmp(s,"ULXMAP") == 0)
-                    map->ulxmap = f;
-                else if(strcmp(s,"ULYMAP") == 0)
-                    map->ulymap = f;
-                else if(strcmp(s,"XDIM") == 0)
-                    map->xdim = f;
-                else if(strcmp(s,"YDIM") == 0)
-                    map->ydim = f;
-            }
-        }
-        free(s);
-        if(fclose(fp))
-        {
-            printf("Error closing file.");
-            return;
-        }
-    }
-    else
-    {
-        printf("Unable to open file %s\n",filename);
-        return;
-    }
-    return;
-}
-
-/**
- * Read data from a single gtopo30 map
- * -- NOTE: Is there a better way to switch byte order? --
- * @param[in] *data - pointer to data array
- * @param[in] *map - pointer to struct containing map information
- * @param[in] *filename - name of file to be read
- */
-void readMap(short int *data, map_info *map, char *filename)
-{
-    FILE *fp = fopen(strcat(filename,".DEM"),"rb");
-    int nrows = map->nrows, ncols = map->ncols, i;
-    short temp;
-    if(fp != NULL)
-    {
-        fread(data,sizeof(short),nrows*ncols,fp);
-        /* Change byte order */
-        for(i=0; i<ncols*nrows; i++)
-        {
-            temp = ntohs(data[i]);
-            data[i] = temp;
-        }
-    }
-    else
-    {
-        printf("Unable to open file %s\n",filename);
-        return;
-    }
-    fclose(fp);
-}
-
-/**
- * Read data from two gtopo30 maps
- * -- NOTE: Is there a better way to switch byte order? --
- * @param[in] *data - pointer to data array
- * @param[in] *map1 - pointer to struct containing map information
- * @param[in] *map2 - pointer to struct containing map information
- * @param[in] *filename1 - name of file to be read
- * @param[in] *filename2 - name of file to be read
- */
-void readMap2(short int *data, map_info *map1, map_info *map2, \
-              char *filename1, char *filename2)
-{
-    FILE *fp1 = fopen(strcat(filename1,".DEM"),"rb");
-    FILE *fp2 = fopen(strcat(filename2,".DEM"),"rb");
-    int i, j;
-    int nrows1 = map1->nrows, ncols1 = map1->ncols, nbytes1 = (map1->nbits)/8;
-    int nrows2 = map2->nrows, ncols2 = map2->ncols, nbytes2 = (map2->nbits)/8;
-    short int temp;
-    short int *data1 = malloc(nrows1*ncols2*nbytes1);
-    short int *data2 = malloc(nrows1*ncols2*nbytes2);
-    if(fp1 != NULL)
-    {
-        fread(data1,sizeof(short),nrows1*ncols1,fp1);
-        /* Change byte order */
-        for(j=0; j<nrows1; j++)
-        {
-            for(i=0; i<ncols1; i++)
-            {
-                temp = ntohs(data1[j*ncols1+i]);
-                data[j*(ncols1+ncols2)+i] = temp;
-            }
-        }
-    }
-    else
-    {
-        printf("Unable to open file %s\n",filename1);
-        return;
-    }
-    free(data1);
-    fclose(fp1);
-    if(fp2 != NULL)
-    {
-        fread(data2,sizeof(short),nrows2*ncols2,fp2);
-        /* Change byte order */
-        for(j=0; j<nrows2; j++)
-        {
-            for(i=0; i<ncols1; i++)
-            {
-                temp = ntohs(data2[j*ncols1+i]);
-                data[j*(ncols1+ncols2)+ncols1+i] = temp;
-            }
-        }
-    }
-    else
-    {
-        printf("Unable to open file %s\n",filename2);
-        return;
-    }
-    free(data2);
-    fclose(fp2);
-}
-
-/**
- * Compute amount of blockage, assuming Gaussian main lobe
- * Blockage is defined by \int_{-elLim}^elMax exp(-(x-b)^2/c) dx / bb_tot
- * where bb_tot is defined through \int_{-elLim}^elLim exp(-x^2/c) dx=1,
- * elLim = sqrt(-c*log(10^(dBlim/10)))), and c=-(beamwidth/2)^2/log(.5);
- * @param[in] *bb - pointer to beam blockage array
- * @param[in] *zi - pointer to interpolated topography array
- * @param[in] lat - latitude of radar in degrees
- * @param[in] height - altitude of radar above sea level in meters
- * @param[in] *range - pointer to radar range in meters
- * @param[in] el - lowest elevation angle of radar in degrees
- * @param[in] beamwidth - beamwidth of radar beam
- * @param[in] dBlim - limit of Gaussian approximation of main lobe
- * @param[in] ri - number of elements in array range
- * @param[in] ai - number of elements in azimuth
- */
-void getBlockage(double *bb, double *zi, double lat, double height, \
-                   double *range, double el, double beamwidth, \
-                   double dBlim, int ri, int ai)
-{
-    /* Declaration of variables */
-    double dndh = -3.9e-8, RE, R, c, elLim, bb_tot;
-    double *elBlock = malloc(sizeof(double)*ri*ai);
-    double *phi = malloc(sizeof(double)*ri);
-    int i, j;
-    
-    RE = getEarthRadius(lat);
-    R = 1./((1./RE) + dndh);
-    for(i=0; i<ai ; i++)
-    {
-        /* Get angle phi from height equation */
-        for(j=0; j<ri ; j++)
-        {
-            phi[j] = rad2deg( asin( ( pow(zi[i*ri+j]+R,2) - pow(range[j],2) - \
-                           pow(R+height,2) ) / (2*range[j]*(R+height)) ) );
-        }
-        /* Find cumulative maximum */
-        cummax(phi, ri);
-        /* Write angles for all area */
-        for(j=0; j<ri ; j++)
-        {
-            elBlock[i*ri+j] = phi[j];
-        }
-    }
-    
-    /* Width of Gaussian */
-    c = -pow(beamwidth/2,2) / log(.5);
-    
-    /* Elevation limits */
-    elLim = sqrt( -c*log( pow(10.,(dBlim/10.)) ) );
-    
-    /* Find total blockage within -elLim to +elLim */
-    bb_tot = sqrt(M_PI*c) * erf(elLim/sqrt(c));
-    
-    for(i=0; i<ai; i++)
-    {
-        for(j=0; j<ri; j++)
-        {
-            /* Check if blockage is seen, otherwise -9999 */
-            if(elBlock[i*ri+j]<el-elLim)
-                elBlock[i*ri+j] = -9999;
-            
-            /* Set maximum for blockage elevation */
-            if(elBlock[i*ri+j]>el+elLim)
-                elBlock[i*ri+j] = el+elLim;
-            
-            bb[i*ri+j] = -1./2. * sqrt(M_PI*c) * ( erf( (el-elBlock[i*ri+j]) / \
-                        sqrt(c) ) - erf( (el-(el-elLim))/sqrt(c) ) ) / bb_tot;
-        }
-    }
-    free(phi);
-    free(elBlock);
-    return;
-}
-
-
diff --git a/lib/beamb_map.h b/lib/beamb_map.h
deleted file mode 100644 (file)
index 2308495..0000000
+++ /dev/null
@@ -1,15 +0,0 @@
-/* Various functions to read map. Used by beamb.c */
-
-int getTopo(double lat, double lon, double max_range);
-
-void readHdr(map_info *map, char *filename);
-
-void readMap(short int *data, map_info *map, char *filename);
-
-void readMap2(short int *data, map_info *map1, map_info *map2, \
-              char *filename1, char *filename2);
-
-void getBlockage(double *bb, double *zi, double lat, double height, \
-                 double *range, double el, double beamwidth, double dBlim, \
-                 int ri, int ai);
-