"""Sensitivity checks for the registered solar relative-release approximation.

Frozen daily JPL states are interpolated with position/velocity cubic Hermite.
This is an independent point-planet sensitivity calculation, not a JPL dust fit.
"""
import math
from comet_orbits import MU, AU_METRES, DAY_SECONDS, add, subtract, scale, dot, magnitude, kepler, position_angle

def interpolated_state(states, time):
    index = max(0, min(len(states)-2, int(math.floor(time-states[0]['jdTdb']))))
    first, second = states[index:index+2]
    interval = second['jdTdb']-first['jdTdb']
    fraction = max(0, min(1, (time-first['jdTdb'])/interval))
    t = fraction
    h00, h10, h01, h11 = 2*t**3-3*t*t+1, t**3-2*t*t+t, -2*t**3+3*t*t, t**3-t*t
    position = [h00*a+h10*interval*b+h01*c+h11*interval*d for a,b,c,d in zip(first['positionAu'],first['velocityAuPerDay'],second['positionAu'],second['velocityAuPerDay'])]
    velocity = [((6*t*t-6*t)*a+(3*t*t-4*t+1)*interval*b+(-6*t*t+6*t)*c+(3*t*t-2*t)*interval*d)/interval for a,b,c,d in zip(first['positionAu'],first['velocityAuPerDay'],second['positionAu'],second['velocityAuPerDay'])]
    return position, velocity

def inverse_cube_field(vector, gravity):
    return scale(vector, -gravity/magnitude(vector)**3)

def cross(first, second):
    return [first[1]*second[2]-first[2]*second[1],first[2]*second[0]-first[0]*second[2],first[0]*second[1]-first[1]*second[0]]

def nuclear_outgassing(position, velocity, parameters):
    radius = magnitude(position)
    radial = scale(position, 1/radius)
    normal = cross(position, velocity)
    normal = scale(normal, 1/magnitude(normal))
    transverse = cross(normal, radial)
    ratio = radius/parameters['R0Au']
    activity = parameters['ALN']*ratio**(-parameters['NM'])*(1+ratio**parameters['NN'])**(-parameters['NK'])
    return scale(add(scale(radial, parameters['A1AuPerDay2']),scale(transverse,parameters['A2AuPerDay2'])),activity)

def relative_orbit(snapshot, age, beta, include_planets=True, include_outgassing=True, step=0.5):
    final_time = snapshot['stopJdTdb']
    start = final_time-age
    state = [0.0]*6
    steps = math.ceil(age/step)
    step = age/steps
    planet_gravity = [row['gravitationalParameterKm3PerS2']*1000**3*DAY_SECONDS**2/AU_METRES**3 for row in snapshot['planets']]
    def derivative(values,time):
        comet, velocity = interpolated_state(snapshot['comet'],time)
        dust = add(comet,values[:3])
        acceleration = subtract(inverse_cube_field(dust,MU*(1-beta)),inverse_cube_field(comet,MU))
        if include_planets:
            for planet,gravity in zip(snapshot['planets'],planet_gravity):
                centre,_ = interpolated_state(planet['states'],time)
                acceleration = add(acceleration,subtract(inverse_cube_field(subtract(dust,centre),gravity),inverse_cube_field(subtract(comet,centre),gravity)))
        if include_outgassing:
            acceleration = subtract(acceleration,nuclear_outgassing(comet,velocity,snapshot['nucleusOutgassing']))
        return values[3:]+acceleration
    for index in range(steps):
        time = start+index*step
        first = derivative(state,time)
        second = derivative(add(state,scale(first,step/2)),time+step/2)
        third = derivative(add(state,scale(second,step/2)),time+step/2)
        fourth = derivative(add(state,scale(third,step)),time+step)
        state = [v+step*(a+2*b+2*c+d)/6 for v,a,b,c,d in zip(state,first,second,third,fourth)]
    return state[:3]

def check_relative_model(snapshot, ephemeris):
    ra = ephemeris['observerAstrometry']['raDegrees']
    dec = ephemeris['observerAstrometry']['decDegrees']
    records = []
    for age in [30,60,90,120]:
        released = min(ephemeris['releaseStates'],key=lambda row:abs(row['ageDays']-age))
        parent,_ = kepler(released['positionAu'],released['velocityAuPerDay'],age)
        for beta in [0.005,0.03,0.1,0.3]:
            dust,_ = kepler(released['positionAu'],released['velocityAuPerDay'],age,MU*(1-beta))
            conic = subtract(dust,parent)
            corrected = relative_orbit(snapshot,age,beta)
            solar_reference = relative_orbit(snapshot,age,beta,False,False)
            planets_without_outgassing = relative_orbit(snapshot,age,beta,True,False)
            difference = magnitude(subtract(corrected,conic))*AU_METRES/1000
            records.append({'ageDays':age,'beta':beta,'conicPositionAngleDegrees':position_angle(conic,ra,dec),
                'perturbedPositionAngleDegrees':position_angle(corrected,ra,dec),
                'totalPositionDifferenceKm':difference,'fractionOfConicSeparation':difference/magnitude(conic)/(AU_METRES/1000),
                'referencePathOnlyDifferenceKm':magnitude(subtract(solar_reference,conic))*AU_METRES/1000,
                'planetAndOutgassingDifferenceKm':magnitude(subtract(corrected,solar_reference))*AU_METRES/1000,
                'unsharedNucleusOutgassingEffectKm':magnitude(subtract(corrected,planets_without_outgassing))*AU_METRES/1000})
    shortest_radius = min(magnitude(row['positionAu']) for row in snapshot['comet'])
    p = snapshot['nucleusOutgassing']
    ratio = shortest_radius/p['R0Au']
    activity = p['ALN']*ratio**(-p['NM'])*(1+ratio**p['NN'])**(-p['NK'])
    bound = 0.5*activity*math.hypot(p['A1AuPerDay2'],p['A2AuPerDay2'])*120**2*AU_METRES/1000
    return {'model':'Solar relative conic releases compared with daily-JPL-path relative RK4, eight point planetary systems and unshared nucleus A1/A2',
            'stepDays':0.5,'records':records,'nucleusMaximumAccelerationFreeDisplacementEstimate120DaysKm':bound,
            'largestPositionAngleDifferenceDegrees':max(abs((row['perturbedPositionAngleDegrees']-row['conicPositionAngleDegrees']+180)%360-180) for row in records),
            'omitted':'No full SPICE/JPL dust fit, planetary multipoles,16 asteroid perturbers or relativistic terms. This sampled sensitivity result is not a rigorous all-parameter error bound.'}