#!/usr/bin/env python
import numpy as np
from   scipy.optimize    import minimize 
#######################################################
def distance(A,B):
        'distance between two points A and B'
        A, B = np.array(A), np.array(B)
        return np.linalg.norm(A-B)
#######################################################

def psi_interp1D_edit(r,Z,E,F):
	
	r_min, r_max = min(r), max(r) + 0.5
	Z_min, Z_max = min(Z), max(Z) - 0.1

	MEP =  [ ] 
	for phi in np.arange(-19., 90., 1.):
#	for phi in [ 84. ]:

		alpha = np.radians(phi)				# interpolation angle
		m = np.tan(alpha)				# slope
		q = Z_max - m * r_max 				# intercept

	# Check boundaries according to phi value ----------------------------------------------
		if -20. < phi <= 45.:
			r_initial       = r_min                 # r and Z range to consider
			r_final         = r_max                 # for the 1D interpolations
			Z_initial       = m*r_min+q		#
			Z_final         = Z_max                 #
		elif 45. < phi < 90.:
			r_initial	= float(Z_min - q) / m
			r_final		= r_max
			Z_initial	= Z_min
			Z_final		= Z_max
		else:
			print( "you should only consider 0 < phi < 90.\nexiting..." )
			quit()
	# --------------------------------------------------------------------------------------

		lenght = distance([r_final, Z_final], [r_initial, Z_initial])   # lenght of the segment in the plot [Ang]
	
		pos_OnTheLine, E_OnTheLine 	= [ ], [ ]	# positions & energy values along the deifined line
		r_OnTheLine, Z_OnTheLine 	= [ ], [ ]	# r,Z points used to evaluate the potential on the line
		position = 0.					# starting point on the line 

		def potential_OnTheLine(percent):			# this function return the value
			x = r_initial+ percent*lenght*np.cos(alpha)	# of the interpolated potential 
			y = y = m * x + q				# on the selected line (known m and q)
			return F(x,y)					# for the percentual position 
									# 0 == beginning, 1 == end	
		minimum = minimize(potential_OnTheLine, 0.15, method = 'Nelder-Mead')
		
		opt_pos = minimum.x					# percentual position of the minimum
		opt_r	= r_initial+opt_pos*lenght*np.cos(alpha)	# minimum r 
		opt_Z 	= m * opt_r + q					# minimum Z
                                 # r mep   # Z mep   # E mep
		MEP.append(     [ opt_r,    opt_Z,   minimum.fun ]      )

	MEP = np.array(MEP)
	return MEP  # return the Minimum energy path