// =====
// Copyright 2007 John Van Vliet, all rights reserved
// Here's my simple license, soon to be replaced by a proper license.
// What you can do:
// - read the code
// - play with the code
// - run the code
// - modify the code
// What you can't do:
// - re-publish the original code or any modified versions
// =====using System;
using System.Collections.Generic;
using System.Text;
namespace OpenTemp
{
/// <summary>
/// Definition of a single station and its weight in determining the temperature for this cell.
/// </summary>
struct CellStation
{
internal float weight;
internal Station station;
}
/// <summary>
/// Definition of a single cell on the earth's surface
/// </summary>
class Cell
{
private const double _earthRadius = 6373f; // approximate radius of earth
private float _latitude = 0.0f;
private float _longitude = 0.0f;
private float _area = 0.0f;
private List<CellStation> _stations = new List<CellStation>();
/// <summary>
/// Constructor
/// </summary>
/// <param name=”latitude”>Latitude of the centre of this cell (in degrees North)</param>
/// <param name=”longitude”>Longitutde of the centre of this cell (in degrees East)</param>
/// <param name=”height”>Height of this cell (in degrees NS)</param>
/// <param name=”width”>Width of this cell (in degrees EW)</param>
internal Cell(float latitude, float longitude, float width, float height)
{
_latitude = latitude;
_longitude = longitude;
// TODO: Improve the area calculation for a cell
// The current calculation works well for a small cell, but will break down for larger cells (> ~1000km on a side)
float sinLatitude = (float)Math.Sin(latitude * Math.PI / 180.0);
float degLength = (float)(2.0 * Math.PI * _earthRadius / 360.0); // Length of a degree at the equator
_area = (width * sinLatitude * degLength) * (height * degLength);
}
/// <summary>
/// Get the latitude (in degrees) of the centre of this cell.
/// </summary>
internal float Latitude
{
get { return _latitude; }
}
/// <summary>
/// Get the longitude (in degrees) of the centre of this cell.
/// </summary>
internal float Longitude
{
get { return _longitude; }
}
/// <summary>
/// Get the area of this cell.
/// </summary>
/// <remarks>The area of the cell is calculated from the width and height. Currently the area calculation assumes
/// small values for width and height. That is, it assumes that the cell is square in cartesian coordinates.
/// The error for a 1000km cell is approximately 0.2%, increasing to 0.8% for a 2000km cell.</remarks>
internal float Area
{
get { return _area; }
}
/// <summary>
/// Get the list of stations that affect the temperature on this cell
/// </summary>
/// <remarks>In the SortedList, the Key is the station weight (0 to 1) and the Value is Station itself.</remarks>
internal List<CellStation> Stations
{
get { return _stations; }
}
/// <summary>
/// Define a list of all stations that may affect the average for this cell.
/// </summary>
/// <param name=”stations”>List of stations to check</param>
/// <param name=”averagingRadius”>Radius of stations to include in calculation</param>
/// <returns>Number of stations that affect the average for this cell</returns>
internal int DefineStations(ICollection<Station> allStations, float averagingRadius)
{
averagingRadius = (float)Math.Abs(averagingRadius);
_stations.Clear();
foreach (Station station in allStations)
{
float stationDistance = CalcDistance(station.Latitude, station.Longitude);
if (stationDistance <= averagingRadius)
{
CellStation cs;
cs.weight = (averagingRadius - stationDistance) / averagingRadius;
cs.station = station;
_stations.Add(cs);
}
}
return _stations.Count;
}
/// <summary>
/// Calculate the distance from this station to a point at the given latitude and longitude.
/// </summary>
/// <param name=”latitude”></param>
/// <param name=”longitude”></param>
/// <returns>Distance from this cell to the given coordinates.</returns>
/// <remarks>The distance calculation uses a straight line between this cell and the given points.
/// It does not follow the curvature of the earth. This simplification results in errors of
/// 0.1% at 1000km, 0.4% at 2000km, and 1.6% at 4000km (all under-estimations).
/// </remarks>
internal float CalcDistance(float latitude, float longitude)
{
double[] posThis = new double[3];
double[] posPoint = new double[3];
posThis[0] = _earthRadius * Math.Sin(_longitude * Math.PI / 180.0) * Math.Cos(_latitude * Math.PI / 180.0);
posThis[2] = _earthRadius * Math.Sin(_latitude * Math.PI / 180.0);
posThis[1] = Math.Sqrt((_earthRadius * _earthRadius) - ((posThis[0] * posThis[0]) + (posThis[2] * posThis[2])));
posPoint[0] = _earthRadius * Math.Sin(longitude * Math.PI / 180.0) * Math.Cos(latitude * Math.PI / 180.0);
posPoint[2] = _earthRadius * Math.Sin(latitude * Math.PI / 180.0);
posPoint[1] = Math.Sqrt((_earthRadius * _earthRadius) - ((posPoint[0] * posPoint[0]) + (posPoint[2] * posPoint[2])));
double dx = posThis[0] - posPoint[0];
double dy = posThis[1] - posPoint[1];
double dz = posThis[2] - posPoint[2];
return (float)Math.Sqrt((dx * dx) + (dy * dy) + (dz * dz));
}
}
}