# Conversion of the (median) location error ellipse for the best location
# of any lightning stroke into the probability that the stroke was within
# any radius of any object of interest. Implementation of the method pre-
# sented in the NASA report [1].
#
# [1] Huddleston LL, Roeder WP, Merceret FJ. A Probabilistic, Facility-Centric
#     Approach to Lightning Strike Location. tech. rep., National Aeronautics
#     and Space Administration; Kennedy Space Center, FL 32899-0001: 2012.
#     NASA/TM-20 12-216308.

# library imports
import numpy as np
from scipy import integrate

# Great circle (orthodromic) distance between points on the Earth's surface
def distance(lon_wt, lat_wt, lon_rad, lat_rad):
    """
    :param lon_wt: longitude of the wind turbine location
    :param lat_wt: latitude of the wind turbine location
    :param lon_rad: longitude of the lightning strike location
    :param lat_rad: latitude of the lightning strike location
    :return: orthodromic distance in meters
    """
    # Convert longitude and latitude information to radians
    lon_wt = lon_wt*np.pi/180.
    lat_wt = lat_wt*np.pi/180.
    lon_rad = lon_rad*np.pi/180.
    lat_rad = lat_rad*np.pi/180.
    # Compute a great circle distance between two points on the globe
    distance = 2.*np.arcsin(np.sqrt((np.sin((lat_wt-lat_rad)/2.))**2 +
                                    np.cos(lat_wt)*np.cos(lat_rad) *
                                    (np.sin((lon_wt-lon_rad)/2.))**2))
    # Convert distance to nautical miles
    distance_nm = ((180.*60.)/np.pi) * distance
    # convert distance to meters
    distance_m = distance_nm * 1852.
    return distance_m

def coordinate_rotation(target_lon, target_lat, strike_lon, strike_lat,
                        alpha, distance, semi_major, semi_minor, proba):
    # convert degrees to radians
    target_lon = target_lon*(np.pi/180.)
    target_lat = target_lat*(np.pi/180.)
    strike_lon = strike_lon*(np.pi/180.)
    strike_lat = strike_lat*(np.pi/180.)
    alpha = alpha*(np.pi/180.)
    # coordinate rotation
    X = (target_lon - strike_lon)*np.cos(strike_lat)
    Z = target_lat - strike_lat
    theta = alpha - (np.pi/2. - np.arctan2(Z, X))
    Xc = distance*np.cos(theta)
    Zc = distance*np.sin(theta)
    k = np.sqrt(-2.*np.log(1. - proba))
    sigma_x = semi_major/k
    sigma_z = semi_minor/k
    return Xc, Zc, sigma_x, sigma_z

def kernel_function(x, mu_k, mu_h, sigma_k, sigma_h, w):
    z1 = (np.sqrt(w-x**2) - mu_k)/(np.sqrt(2.)*sigma_k)
    z2 = (np.sqrt(w-x**2) + mu_k)/(np.sqrt(2.)*sigma_k)
    erf_z1 = special.erf(z1)
    erf_z2 = special.erf(z2)
    exp_1 = np.exp(-(x - mu_h)**2/(2.*sigma_h**2))
    exp_2 = np.exp(-(x + mu_h)**2/(2.*sigma_h**2))
    function = (exp_1 + exp_2)*(erf_z1 + erf_z2)
    return function

# EXAMPLE:
# lightning strike location
strike_lat = 28.6069
strike_lon = -80.6087
# confidence ellipse (50%)
semi_major = 0.6  # km
semi_minor = 0.4  # km
heading = 82.  # heading from true north
# point of interest (e.g. WT location)
poi_lat = 28.60827
poi_lon = -80.6041
# radius around point of interest
radius = 0.45*1.852  # nautical miles to km

# COMPUTATION:
# orthodromic distance
dist = distance(poi_lon, poi_lat, strike_lon, strike_lat)
dist = dist*1e-3  # m => km
# coordinate system rotation
proba = 0.5  # confidence ellipse probability
mu_k, mu_h, sigma_k, sigma_h = coordinate_rotation(poi_lon, poi_lat, strike_lon, strike_lat,
                                                   heading, dist, semi_major, semi_minor, proba)
w = radius**2
arguments = (mu_k, mu_h, sigma_k, sigma_h, w)
# numerical integration (quadrature)
integral = integrate.quad(kernel_function, 0., np.sqrt(w), args=arguments)
# probability of strike within a circle
konstanta = 1./(2.*np.sqrt(2.*np.pi)*sigma_h)
P = konstanta*integral[0]
# probability that a strike is within the radius around the point of interest
print(P)