/* proj.c */
/*
* Vis5D system for visualizing five dimensional gridded data sets.
* Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul,
* Dave Santek, and Andre Battaiola.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* As a special exception to the terms of the GNU General Public
* License, you are permitted to link Vis5D with (and distribute the
* resulting source and executables) the LUI library (copyright by
* Stellar Computer Inc. and licensed for distribution with Vis5D),
* the McIDAS library, and/or the NetCDF library, where those
* libraries are governed by the terms of their own licenses.
*
* This program 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "../config.h"
/*
* Map projection / coordinate transformation
*/
/*
* Map projection parameters:
*
* IF Projection==PROJ_GENERIC THEN
* projargs[0] = northbound
* projargs[1] = westbound
* projargs[2] = row increment
* projargs[3] = column increment
* ELSE IF Projection==PROJ_LINEAR THEN
* projargs[0] = northbound in degrees
* projargs[1] = westbound in deg.
* projargs[2] = row increment in deg.
* projargs[3] = column increment in deg.
* ELSE IF Projection==PROJ_LAMBERT THEN
* projargs[0] = Lat1
* projargs[1] = Lat2
* projargs[2] = PoleRow
* projargs[3] = PoleColumn
* projargs[4] = Central Long.
* projargs[5] = Column Increment in km
* ELSE IF Projection==PROJ_STEREO THEN
* projargs[0] = Central Latitude
* projargs[1] = Central Longitude
* projargs[2] = Central Row
* projargs[3] = Central Column
* projargs[4] = Column increment in km
* ELSE IF Projection==PROJ_ROTATED THEN
* projargs[0] = Latitude on rotated globe of grid row 0
* projargs[1] = Longitude on rotated globe of grid column 0
* projargs[2] = degs of latitude on rotated globe between grid rows
* projargs[3] = degs of longitude on rotated globe between grid cols
* projargs[4] = Earth latitude of (0, 0) on rotated globe
* projargs[5] = Earth longitude of (0, 0) on rotated globe
* projargs[6] = Clockwise rotation of rotated globe in degrees
* ELSE IF Projection==PROJ_MERCATOR THEN
* projargs[0] = Central Latitude
* projargs[1] = Central Longitude
* projargs[2] = Row Increment in Kilometers
* projargs[3] = Column Increment in Kilometers
* ELSE IF Projection==PROJ_EPA THEN
* projargs[ (R*nc+C) * 2 + 0 ] = latitude of grid point (R,C)
* projargs[ (R*nc+C) * 2 + 1 ] = longitude of grid point (R,C)
* ENDIF
*/
/*
* Vertical coordinate system parameters:
*
* IF VertSystem==VERT_GENERIC THEN
* vertargs[0] = Bottom bound
* vertargs[1] = level increment
* ELSE IF VertSystem==VERT_EQUAL_KM THEN
* vertargs[0] = bottom bound in km
* vertargs[1] = level increment in km
* ELSE IF VertSystem==VERT_UNEQUAL_KM THEN
* vertargs[n] = height of grid level n. n is in [0..Nl-1]
* ELSE IF VertSystem==VERT_UNEQUAL_MB THEN
* vertargs[n] = height of grid level n. n is in [0..Nl-1]
* ELSE IF VertSystem==VERT_EPA THEN
* vertargs[n] = height of grid level n. n is in [0..Nl-1]
* ENDIF
*/
#include <math.h>
#include <stdio.h>
#include "grid_i.h"
#include "proj_i.h"
#include "v5d.h"
#ifndef M_PI
# define M_PI 3.1415926
#endif
#define DEG2RAD (M_PI/180.0)
#define RAD2DEG (180.0/M_PI)
#define RADIUS 6371.23
/* Pete rotated sphere to Earth */
static void pandg_back( float *lat, float *lon, float a, float b, float r )
{
float pr, gr, pm, gm;
/* NOTE - longitude sign switches - b too! */
pr = DEG2RAD * *lat;
gr = -DEG2RAD * *lon;
pm = asin( cos(pr) * cos (gr) );
gm = atan2(cos(pr) * sin (gr), -sin(pr) );
*lat = RAD2DEG * asin( sin(a) * sin(pm) - cos(a) * cos(pm) * cos (gm - r) );
*lon = -RAD2DEG * (-b + atan2(cos(pm) * sin (gm - r),
sin(a) * cos(pm) * cos (gm - r) + cos(a) * sin(pm)));
}
/* Earth to Pete rotated sphere */
static void pandg_for( float *lat, float *lon, float a, float b, float r )
{
float p1, g1, p, g;
/* NOTE - longitude sign switches - b too! */
p1 = DEG2RAD * *lat;
g1 = -DEG2RAD * *lon;
p = asin( sin(a) * sin(p1) + cos(a) * cos(p1) * cos (g1 + b) );
g = r + atan2(cos(p1) * sin (g1 + b),
sin(a) * cos(p1) * cos (g1 + b) - cos(a) * sin(p1) );
*lat = RAD2DEG * asin( -cos(p) * cos (g) );
*lon = -RAD2DEG * atan2(cos(p) * sin (g), sin(p) );
}
/**********************************************************************/
/***** Map Projection *****/
/**********************************************************************/
/*
* Convert a (row,column) grid coordinate to a (lat,lon) coordinate
* according to the projection info.
* Return: 1 = ok, 0 = error (out of bounds)
*/
int rowcol_to_latlon_i( float row, float col, float *lat, float *lon,
struct projection *proj )
{
switch (proj->Kind) {
case PROJ_GENERIC: /* Generic rectilinear */
case PROJ_LINEAR: /* lat/lon rectilinear */
case PROJ_CYLINDRICAL:
case PROJ_SPHERICAL:
/* FALL-THROUGH */
#define NorthBound proj->Args[0]
#define WestBound proj->Args[1]
#define RowInc proj->Args[2]
#define ColInc proj->Args[3]
*lat = NorthBound - row * RowInc;
*lon = WestBound - col * ColInc;
#undef NorthBound
#undef WestBound
#undef RowInc
#undef ColInc
break;
case PROJ_LAMBERT: /* Lambert conformal */
#define PoleRow proj->Args[2]
#define PoleCol proj->Args[3]
#define CentralLon proj->Args[4]
#define Hemisphere proj->AuxArgs[0]
#define ConeFactor proj->AuxArgs[1]
#define Cone proj->AuxArgs[2]
{
float xldif, xedif, xrlon, radius;
xldif = Hemisphere * (row-PoleRow) / ConeFactor;
xedif = (PoleCol-col) / ConeFactor;
if (xldif==0.0 && xedif==0.0)
xrlon = 0.0;
else
xrlon = atan2( xedif, xldif );
*lon = xrlon / Cone * RAD2DEG + CentralLon;
if (*lon > 180.0)
*lon -= 360.0;
radius = sqrt( xldif*xldif + xedif*xedif );
if (radius < 0.0001)
*lat = 90.0 * Hemisphere; /* +/-90 */
else
*lat = Hemisphere
* (90.0 - 2.0*atan(exp(log(radius)/Cone))*RAD2DEG);
}
#undef PoleRow
#undef PoleCol
#undef CentralLon
#undef Hemisphere
#undef ConeFactor
#undef Cone
break;
case PROJ_MERCATOR: /* Mercator */
#define CentralLat proj->Args[0]
#define CentralLon proj->Args[1]
#define RowIncKm proj->Args[2]
#define ColIncKm proj->Args[3]
{
float YC, alpha, ic, jc;
YC = RADIUS * log((1.0 + sin(DEG2RAD*CentralLat))/cos(DEG2RAD*CentralLat));
ic = (proj->Nr-1)/2.0;
jc = (proj->Nc-1)/2.0;
alpha = ( (ic-row) * RowIncKm + YC) / RADIUS;
*lat = 2 * RAD2DEG * atan( exp(alpha) ) - 90.0;
*lon = CentralLon - RAD2DEG * (col-jc) * ColIncKm / RADIUS;
}
#undef CentralLat
#undef CentralLon
#undef RowIncKm
#undef ColIncKm
break;
case PROJ_STEREO: /* Polar Stereographic */
#define CentralLat proj->Args[0]
#define CentralLon proj->Args[1]
#define CentralRow proj->Args[2]
#define CentralCol proj->Args[3]
#define ColInc proj->Args[4]
{
float CosCentralLat = cos( CentralLat * DEG2RAD );
float SinCentralLat = sin( CentralLat * DEG2RAD );
float Scale = (2.0 * RADIUS / ColInc);
float InvScale = 1.0 / Scale;
float xrow, xcol, rho, c, cc, sc;
xrow = CentralRow - row;
xcol = CentralCol - col;
rho = xrow*xrow + xcol*xcol;
if (rho<1.0e-20) {
*lat = CentralLat;
*lon = CentralLon;
}
else {
rho = sqrt( rho );
c = 2.0 * atan( rho * InvScale);
cc = cos(c);
sc = sin(c);
*lat = RAD2DEG
* asin( cc*SinCentralLat + xrow*sc*CosCentralLat / rho );
*lon = CentralLon + RAD2DEG * atan2( xcol * sc,
(rho * CosCentralLat * cc - xrow * SinCentralLat * sc) );
if (*lon < -180.0) *lon += 360.0;
else if (*lon > 180.0) *lon -= 360.0;
}
}
#undef CentralLat
#undef CentralLon
#undef CentralRow
#undef CentralCol
#undef ColInc
break;
case PROJ_ROTATED: /* Rotated */
#define NorthBound proj->Args[0]
#define SouthBound (proj->Args[0] - proj->Args[2]*(proj->Nr-1))
#define WestBound proj->Args[1]
#define EastBound (proj->Args[1] - proj->Args[3]*(proj->Nc-1))
#define CentralLat proj->Args[4]
#define CentralLon proj->Args[5]
#define Rotation (proj->Args[6] * DEG2RAD)
*lat = NorthBound - row * (NorthBound-SouthBound) / (float) (proj->Nr-1);
*lon = WestBound - col * (WestBound-EastBound) / (float) (proj->Nc-1);
pandg_back(lat, lon, CentralLat, CentralLon, Rotation);
#undef NorthBound
#undef SouthBound
#undef WestBound
#undef EastBound
#undef CentralLat
#undef CentralLon
#undef Rotation
break;
case PROJ_EPA: /* EPA projection */
#define Latitude(ROW,COL) proj->Args[ ((ROW)*nc+(COL)) * 2 + 0 ]
#define Longitude(ROW,COL) proj->Args[ ((ROW)*nc+(COL)) * 2 + 1 ]
{
/* TODO: interpolate between grid points */
int nr, nc;
nr = proj->Nr;
nc = proj->Nc;
*lat = Latitude( (int) row, (int) col );
*lon = Longitude( (int) row, (int) col );
}
#undef Latitude
#undef Longitude
break;
default:
printf("Error in rowcol_to_latlon_i: bad projection: %d\n",
(int)proj );
break;
}
return 1;
}
/*
* Find the rectangle in the EPA Lambert conformal grid which contains the
* given lat/lon coordinate.
* Input: lat, lon - coordinate of interest.
* proj - the projection
* Output: n, m - row, column of upper-left corner of rectangle.
* s, t - coefficients for interpolation.
* Return: 1 = found a containing rectangle
* 0 = didn't find a rectangle.
*/
static int find_rectangle( float lat, float lon, int *n, int *m,
float *s, float *t, struct projection *proj )
{
#define Latitude(ROW,COL) proj->Args[ ((ROW)*nc+(COL)) * 2 + 0 ]
#define Longitude(ROW,COL) proj->Args[ ((ROW)*nc+(COL)) * 2 + 1 ]
int i, j;
float px, py, dx, dy;
int nr, nc;
/* First two elements of projargs indicate no. of rows and columns */
nr = proj->Nr;
nc = proj->Nc;
for (i=0;i<nr-1;i++) {
for (j=0;j<nc-1;j++) {
/*
* Test if (lat,lon) is inside the rectangle by constructing
* edge vectors, computing the cross product, and testing if
* the z component is always positive.
*/
/* top edge */
px = lon - Longitude(i,j);
py = lat - Latitude(i,j);
dx = Longitude(i,j+1) - Longitude(i,j);
dy = Latitude(i,j+1) - Latitude(i,j);
if (dx*py-dy*px<0.0) {
/* outside */
/* printf("above\n");*/
continue;
}
/* right edge */
px = lon - Longitude(i,j+1);
py = lat - Latitude(i,j+1);
dx = Longitude(i+1,j+1) - Longitude(i,j+1);
dy = Latitude(i+1,j+1) - Latitude(i,j+1);
if (dx*py-dy*px<0.0) {
/* outside */
/*printf("right\n");*/
continue;
}
/* bottom edge */
px = lon - Longitude(i+1,j+1);
py = lat - Latitude(i+1,j+1);
dx = Longitude(i+1,j) - Longitude(i+1,j+1);
dy = Latitude(i+1,j) - Latitude(i+1,j+1);
if (dx*py-dy*px<0.0) {
/* outside */
/*printf("below\n");*/
continue;
}
/* left edge */
px = lon - Longitude(i+1,j);
py = lat - Latitude(i+1,j);
dx = Longitude(i,j) - Longitude(i+1,j);
dy = Latitude(i,j) - Latitude(i+1,j);
if (dx*py-dy*px<0.0) {
/* outside */
/* printf("left\n");*/
continue;
}
/* if we get here, all edge tests succeeded, */
/* therefore, (lat,lon) is inside the rectangle */
/* return row, column location */
*n = i;
*m = j;
/* compute interpolation coefficients */
*s = (lat - Latitude(i,j)) / (Latitude(i+1,j) - Latitude(i,j));
*t = (lon - Longitude(i,j)) / (Longitude(i,j+1) - Longitude(i,j));
return 1;
}
}
return 0;
#undef Latitude
#undef Longitude
}
/*
* Convert a lat/lon coordinate to a row/column according to the given
* projection.
* Input: lat, lon - the lat/lon coordinate
* proj - which projection
* projargs - projection arguments
* nr, nc - number of rows, columns in corresponding 3-D grid
* Output: row, col - the computed row and column
* Return: 1 = ok, 0 = error (lat/lon out of domain)
*/
int latlon_to_rowcol_i( float lat, float lon, float *row, float *col,
struct projection *proj )
{
float lat0, lon0;
switch (proj->Kind) {
case PROJ_GENERIC: /* Generic rectilinear */
case PROJ_LINEAR: /* lat/lon rectilinear */
case PROJ_CYLINDRICAL:
case PROJ_SPHERICAL:
#define NorthBound proj->Args[0]
#define WestBound proj->Args[1]
#define RowInc proj->Args[2]
#define ColInc proj->Args[3]
*row = (NorthBound - lat) / RowInc;
*col = (WestBound - lon) / ColInc;
if (*row<0.0 || *row>(proj->Nr-1) || *col<0.0 || *col>(proj->Nc-1)) {
return 0;
}
#undef NorthBound
#undef WestBound
#undef RowInc
#undef ColInc
break;
case PROJ_LAMBERT: /* Lambert conformal */
#define PoleRow proj->Args[2]
#define PoleCol proj->Args[3]
#define CentralLon proj->Args[4]
#define Hemisphere proj->AuxArgs[0]
#define ConeFactor proj->AuxArgs[1]
#define Cone proj->AuxArgs[2]
{
float rlon, rlat, r;
rlon = lon - CentralLon;
rlon = rlon * Cone * DEG2RAD;
if (lat<-85.0) {
/* infinity */
r = 10000.0;
}
else {
rlat = (90.0 - Hemisphere * lat) * DEG2RAD * 0.5;
r = ConeFactor * pow( tan(rlat), Cone );
}
*row = PoleRow + r * cos(rlon);
*col = PoleCol - r * sin(rlon);
if ( *row<0.0 || *row>(proj->Nr-1)
|| *col<0.0 || *col>(proj->Nc-1)) {
return 0;
}
}
#undef PoleRow
#undef PoleCol
#undef CentralLon
#undef Hemisphere
#undef ConeFactor
#undef Cone
break;
case PROJ_MERCATOR: /* Mercator */
#define CentralLat proj->Args[0]
#define CentralLon proj->Args[1]
#define RowIncKm proj->Args[2]
#define ColIncKm proj->Args[3]
{
float X, Y, ic, jc, YC;
ic = (proj->Nr-1)/2.0;
jc = (proj->Nc-1) / 2.0;
YC = RADIUS * log((1.0 + sin(DEG2RAD*CentralLat))/cos(DEG2RAD*CentralLat));
X = RADIUS * (lon-CentralLon) / RAD2DEG;
Y = RADIUS * log( (1+sin(DEG2RAD*lat))/cos(DEG2RAD*lat) );
*row = ic - (Y - YC)/RowIncKm;
*col = jc - X/ColIncKm;
}
#undef CentralLat
#undef CentralLon
#undef RowIncKm
#undef ColIncKm
break;
case PROJ_STEREO: /* Polar Stereographic */
#define CentralLat proj->Args[0]
#define CentralLon proj->Args[1]
#define CentralRow proj->Args[2]
#define CentralCol proj->Args[3]
#define ColInc proj->Args[4]
{
float CosCentralLat = cos( CentralLat * DEG2RAD );
float SinCentralLat = sin( CentralLat * DEG2RAD );
float Scale = (2.0 * RADIUS / ColInc);
float rlat, rlon, clon, clat, k;
rlat = DEG2RAD * lat;
rlon = DEG2RAD * (CentralLon - lon);
clon = cos(rlon);
clat = cos(rlat);
k = Scale
/ (1.0 + SinCentralLat*sin(rlat) + CosCentralLat*clat*clon);
*col = CentralCol + k * clat * sin(rlon);
*row = proj->Nr - CentralRow
- k * (CosCentralLat * sin(rlat) - SinCentralLat * clat * clon);
if ( *row<0.0 || *row>(proj->Nr-1)
|| *col<0.0 || *col>(proj->Nc-1)) {
return 0;
}
}
#undef CentralLat
#undef CentralLon
#undef CentralRow
#undef CentralCol
#undef ColInc
#undef Nr
break;
case PROJ_ROTATED:
#define NorthBound proj->Args[0]
#define SouthBound (proj->Args[0] - proj->Args[2]*(proj->Nr-1))
#define WestBound proj->Args[1]
#define EastBound (proj->Args[1] - proj->Args[3]*(proj->Nc-1))
#define CentralLat proj->Args[4]
#define CentralLon proj->Args[5]
#define Rotation (proj->Args[6] * DEG2RAD)
lat0 = lat;
lon0 = lon;
pandg_for( &lat0, &lon0, CentralLat, CentralLon, Rotation );
*col = (WestBound-lon0)/proj->Args[3];
*row = (NorthBound-lat0) / proj->Args[2];
if ( *row<0.0 || *row>(proj->Nr-1)
|| *col<0.0 || *col>(proj->Nc-1)) {
return 0;
}
#undef NorthBound
#undef SouthBound
#undef WestBound
#undef EastBound
#undef CentralLat
#undef CentralLon
#undef Rotation
break;
case PROJ_EPA: /* EPA projection */
{
int r, c;
float alpha, beta;
if (find_rectangle( lat, lon, &r, &c, &alpha, &beta, proj )) {
*row = (float) r + alpha;
*col = (float) c + beta;
}
else {
return 0;
}
}
break;
default:
printf("Error in latlon_to_rowcol_i: bad projection: %d\n",
(int)proj );
break;
}
return 1;
}
/*
* Return a floating point value indicitive of the resolution of a
* map projection. Low values indicate high-resolution, high values
* indicate low-resolution. Typically, this value is the product of
* the row and column increments.
* Input: proj - the map projection
* Return: resolution value
*/
float proj_resolution( struct projection *proj )
{
float res;
switch (proj->Kind) {
case PROJ_GENERIC: /* Generic rectilinear */
case PROJ_LINEAR: /* lat/lon rectilinear */
case PROJ_CYLINDRICAL:
case PROJ_SPHERICAL:
#define RowInc proj->Args[2]
#define ColInc proj->Args[3]
res = RowInc * ColInc;
#undef RowInc
#undef ColInc
break;
case PROJ_MERCATOR:
#define RowInc proj->Args[2]/111.75
#define ColInc proj->Args[3]/111.75
res = RowInc * ColInc;
#undef RowInc
#undef ColInc
break;
case PROJ_ROTATED:
#define RowInc proj->Args[2]
#define ColInc proj->Args[3]
res = RowInc * ColInc;
#undef RowInc
#undef ColInc
break;
case PROJ_LAMBERT: /* Lambert conformal */
case PROJ_STEREO: /* Polar Stereographic */
case PROJ_EPA: /* EPA projection */
/* measure distance between a few grid points in center of grid */
{
float midrow = proj->Nr / 2.0;
float midcol = proj->Nc / 2.0;
float lat0, lon0, lat1, lon1, dlat, dlon;
rowcol_to_latlon_i( midrow, midcol, &lat0, &lon0, proj );
rowcol_to_latlon_i( midrow+1.0, midcol+1.0, &lat1, &lon1, proj );
dlat = lat1-lat0;
dlon = lon1-lon0;
res = sqrt( dlat*dlat + dlon*dlon );
}
break;
default:
printf("Error in proj_resolution: bad projection: %d\n", (int)proj );
break;
}
/* Return absolute value */
if (res<0.0) {
return -res;
}
else {
return res;
}
}
/**********************************************************************/
/***** Vertical Coordinate System *****/
/**********************************************************************/
/*
* Perform a binary search of the array for the value and return
* the index as a float. Return a -1.0 if the value is out of bounds.
*/
static float binary_search( float value, float array[], int size )
{
int low, high, mid;
float x;
if (value==array[0]) {
return 0.0;
}
else if (value<array[0]) {
return -1.0;
}
else if (value>array[size-1]) {
return -1.0;
}
else if (value==array[size-1]) {
return (float) (size-1);
}
else {
/* do a binary search of array[] for value */
low = 0;
high = size-1;
while (low<=high) {
mid = (low+high)/2;
if (value<array[mid])
high = mid - 1;
else if (value>array[mid])
low = mid + 1;
else
return (float) mid;
}
/* interpolate a value between high and low */
x = (value-array[high]) / (array[low]-array[high]);
return high * (1.0-x) + low * x;
}
}
/*
* Convert a grid level into a height value.
* Input: level - the grid level
* vcs - the vertical coordinate system
* topo_elev - elevation of topography at the implicit lat/lon
* (only used for EPA projection)
* Output: height - the height value
*/
int level_to_height( float level, float *height, struct vcs *vcs,
float topo_elev )
{
switch (vcs->Kind) {
case VERT_GENERIC:
*height = vcs->Args[0] + level * vcs->Args[1];
break;
case VERT_EQUAL_KM:
*height = vcs->Args[0] + level * vcs->Args[1];
break;
case VERT_UNEQUAL_MB: /* Unequally spaced mb */
case VERT_UNEQUAL_KM: /* Unequally spaced km */
{
int ilevel, ilevel_next;
float alpha;
ilevel = ((int) level < vcs->Nl) ? ((int) level) : (vcs->Nl-1);
ilevel_next = (ilevel+1 < vcs->Nl) ? (ilevel+1) : (vcs->Nl-1);
alpha = level - ilevel;
*height = vcs->Args[ilevel] * (1.0-alpha)
+ vcs->Args[ilevel_next] * alpha;
}
break;
case VERT_EPA:
#define Sigma(K) vcs->Args[(K)]
{
float psurf, s, p;
psurf = 1012.5 * exp( -topo_elev / 7.2 );
s = Sigma( (int) level );
p = 100.0 + s * (psurf - 100.0);
*height = -7.2 * log( p / 1012.5 );
}
#undef Sigma
break;
default: /* error */
printf("Error in level_to_height, can't handle vcs kind %d\n",
vcs->Kind);
}
return 1;
}
/*
* Convert a height to a grid level.
* Input: height - the input height value
* vcs - the vertical coordinate system
* topo_elev - elevation of the topograph at the lat/lon of interest
* nl - number of levels in input coordinate system grid
* Output: level - the output grid level
* Return: 1 = ok, 0 = error (height outside of domain)
*/
int height_to_level( float height, float *level, struct vcs *vcs,
float topo_elev )
{
float lev;
switch (vcs->Kind) {
case VERT_GENERIC:
case VERT_EQUAL_KM:
lev = (height - vcs->Args[0]) / vcs->Args[1];
break;
case VERT_UNEQUAL_MB: /* Unequally spaced mb */
case VERT_UNEQUAL_KM: /* Unequally spaced km */
/* lev will be -1 if out of bounds */
lev = binary_search( height, vcs->Args, vcs->Nl );
break;
case VERT_EPA:
#define Sigma(K) vcs->Args[(K)]
{
float hgt[MAXLEVELS];
int ilev;
float psurf, s, p;
/* compute hgt[ilev] */
psurf = 1012.5 * exp( -topo_elev / 7.2 );
for (ilev=0;ilev<vcs->Nl;ilev++) {
s = Sigma( ilev );
p = 100.0 + s * (psurf - 100.0);
hgt[ilev] = -7.2 * log( p / 1012.5 );
}
/* search hgt[] for height */
/* lev will be -1 if out of bounds */
lev = binary_search( height, hgt, vcs->Nl );
#undef Sigma
}
break;
default: /* error */
printf("Error in height_to_level, can't handle vcs kind %d\n",
vcs->Kind);
}
/* boundary checks: */
if ((lev < vcs->LowLev) || (lev > (vcs->Nl-1))) {
return 0;
}
else {
*level = lev;
return 1;
}
}
syntax highlighted by Code2HTML, v. 0.9.1