/*
* EasyWave - A realtime tsunami simulation program with GPU support.
* Copyright (C) 2014 Andrey Babeyko, Johannes Spazier
* GFZ German Research Centre for Geosciences (http://www.gfz-potsdam.de)
*
* Parts of this program (especially the GPU extension) were developed
* within the context of the following publicly funded project:
* - TRIDEC, EU 7th Framework Programme, Grant Agreement 258723
* (http://www.tridec-online.eu)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*/
#include
#include
#include
#include
#include
//=========================================================================
// Constructor
cObsArray::cObsArray()
{
nPos = 0;
nObs = 0;
id = NULL;
lon = NULL;
lat = NULL;
obs = NULL;
}
//=========================================================================
// Destructor
cObsArray::~cObsArray()
{
if( lon ) delete [] lon;
if( lat ) delete [] lat;
if( obs ) {
for( int n=0; n 0 ) {
obs = new double* [nPos];
for( n=0; n 1 ) rad = 1;
c = 2 * asin(rad);
dist = REARTH*c;
return dist;
}
double GeoStrikeOnSphere( double lon1, double lat1, double lon2, double lat2 )
{
const double G2R = 3.14159265358979/180.; // multiplyer to convert from degrees into radians
const double R2G = 180./3.14159265358979; // multiplyer to convert from radians into degrees
const double REARTH = 6378.137; // Earth radius in km along equator
double strike, distRad;
if( (lat1 == lat2) && (lon1 == lon2) ) {
strike = 0.;
}
else if( lon1 == lon2 ) {
if( lat1 > lat2 )
strike = 180.;
else
strike = 0.;
}
else {
distRad = GeoDistOnSphere( lon1,lat1,lon2,lat2 ) / REARTH;
strike = R2G * asin( cos(G2R*lat2)*sin(G2R*(lon2-lon1)) / sin(distRad) );
if( (lat2 > lat1) && (lon2 > lon1) ) {
}
else if( (lat2 < lat1) && (lon2 < lon1) ) {
strike = 180.0 - strike;
}
else if( (lat2 < lat1) && (lon2 > lon1) ) {
strike = 180.0 - strike;
}
else if( (lat2 > lat1) && (lon2 < lon1) ) {
strike += 360.0;
}
}
// if( strike > 180.0 ) strike -= 180.0;
return strike;
}