#Written in Python 3.4.4

#Heave motion calculation script
#Added mass and drag coefficients calculated using DNV RP H103 recommended practice document.

#Module import
import configparser
import numpy as np
import scipy, math
import pandas
import datetime as dt
from pandas import DataFrame, Series
from scipy.interpolate import interp1d as inter
from scipy.optimize import curve_fit as cf
from scipy import constants
from sympy import sin, cos, exp, integrate, simplify, expand_trig, Symbol, symbols
import matplotlib.pyplot as plot

#float range function for use in specifying frequency range
def frange(start, stop, n):
    L = [0.0] * n
    divs = n - 1
    divsinv = 1.0 / divs
    for i in range(n):
        L[i] = divsinv * (start*(divs - i) + stop*i)
    return L		

#function to create wave spectrum for standard sea states (bretschneider modified PM)
def pm(sea_state, freq):
	mean_sig_wave_height = [0.05, 0.05, 0.3, 0.88, 1.88, 3.25, 5, 7.5, 11.5, 14]
	mp_wp = [0, 0, 7.5, 7.5, 8.8, 9.7, 12.4, 15, 16.4, 20]
	try:
		wp = 2 * math.pi / mp_wp[sea_state]		
		return (5/16)* mean_sig_wave_height[sea_state]**2 * wp**4 * (1/freq**5) * exp(-(5/4)*(1/((freq/wp)**4)))
	except IndexError:
		print ('Invalid Sea State. Please choose a sea state between 0 and 9.')
		
def pm_usrdef(Hs, Tp, freq):
	wp = 2 * math.pi / Tp		
	return (5/16)* Hs**2 * wp**4 * (1/freq**5) * exp(-(5/4)*(1/((freq/wp)**4)))
	
		

#declare basic variables

config = configparser.ConfigParser()
config.read('heavecalc_settings.ini')
		
#DDS Geometry
dds_draft = float(config['deep_draft_semi']['draft'])
dds_length = float(config['deep_draft_semi']['length'])
dds_pontoon_width = float(config['deep_draft_semi']['pontoon_width'])
dds_pontoon_height = float(config['deep_draft_semi']['pontoon_height'])
dds_square_leg_width = float(config['deep_draft_semi']['square_column_width'])
dds_pontoon_length = dds_length-(2*dds_square_leg_width)

#EDS Geometry
eds_draft = float(config['deep_draft_semi']['draft'])
eds_length = float(config['deep_draft_semi']['length'])
eds_pontoon_width = float(config['deep_draft_semi']['pontoon_width'])
eds_pontoon_height = float(config['deep_draft_semi']['pontoon_height'])
eds_square_leg_width = float(config['deep_draft_semi']['square_column_width'])
eds_pontoon_length = eds_length-(2*eds_square_leg_width)

heave_plate_separation = float(config['extendable_draft_semi']['heave_plate_separation'])
heave_plate_cutout_width = float(config['extendable_draft_semi']['heave_plate_cutout_width'])
heave_plate_height = float(config['extendable_draft_semi']['heave_plate_height'])
heave_plate_section_width = (eds_length - heave_plate_cutout_width) / 2

extensor_diameter = float(config['extendable_draft_semi']['extensor_diameter'])
extensor_sep = float(config['extendable_draft_semi']['extensor_separation'])
extensor_spacing = float(config['extendable_draft_semi']['extensor_spacing'])

#range values for varying heave plate depth
sep_s = float(config['eds_vary_heave_plate_depth']['start_sep'])
sep_e = float(config['eds_vary_heave_plate_depth']['end_sep'])
sep_divs = int(config['eds_vary_heave_plate_depth']['n_divs'])

#Wave Data
rho = float(config['wave_data']['rho'])
zeta = float(config['wave_data']['damping_ratio'])
sea_state = int(config['wave_data']['sea_state'])
Hs_user = float(config['wave_data']['userdef_Hs'])
Tp_user = float(config['wave_data']['userdef_Tp'])

#Local Constants
g = scipy.constants.value(u'standard acceleration of gravity')

'''
Handling Booleans
'''

dds_do = int(config['deep_draft_semi']['perform_dds_analysis']) == 1
eds_do = int(config['extendable_draft_semi']['perform_eds_analysis']) == 1
var_do = int(config['eds_vary_heave_plate_depth']['perform_varying_eds_analysis']) == 1

pms_do = int(config['wave_data']['usestd']) == 1

#single results
ret_mass = int(config['results']['return_mass']) == 1
ret_addedmass = int(config['results']['return_added_mass']) == 1
ret_wn = int(config['results']['return_natural_frequency']) == 1
ret_tn = int(config['results']['return_natural_period']) == 1

#tabular results
	#format
ret_plot = int(config['results']['return_as_plot']) == 1
ret_x_w = int(config['results']['plot_freq_xaxis']) == 1
ret_xls = int(config['results']['return_as_xls']) == 1
	#series
ret_RAO = int(config['results']['return_RAO']) == 1	
ret_response = int(config['results']['return_response_spectrum']) == 1
ret_spec = int(config['results']['include_wave_spec']) == 1

#characteristic values

ret_sig = int(config['results']['return_sig_heave_amp']) == 1
ret_mean = int(config['results']['return_mean_heave_amp']) == 1
ret_max = int(config['results']['return_max_heave_amp']) == 1

#wave data
ang_freq = Symbol('w')
w_num = Symbol('k')
w_amp = 1
z, x, t = symbols('z, x, t')

w_l = [x for x in frange(float(config['wave_data']['freq_rads_start']), float(config['wave_data']['freq_rads_stop']), int(config['wave_data']['freq_ndiv'])+1)]
period = [2*math.pi/x for x in w_l]
k = [x**2/g for x in w_l]
if pms_do:
	spec = [pm(sea_state, x) for x in w_l]
elif not pms_do:
	spec = [pm_usrdef(Hs_user, Tp_user, x) for x in w_l]	

#2d rectangular ca
ab2dr = [2**32, 10, 5, 2, 1, 0.5, 0.2, 0.1]
ca2dr = [1, 1.14, 1.21, 1.36, 1.51, 1.7, 1.98, 2.23]
try:
	twod_cafunc = inter(ab2dr, ca2dr, kind ='cubic')
except np.linalg.linalg.LinAlgError:	
	twod_cafunc = inter(ab2dr, ca2dr, kind ='quadratic')

#3D rectangular plate added mass coefficients as per DNV RP H103
barect = np.array([1, 1.25, 1.5, 1.59, 2, 2.5, 3, 3.17, 4, 5, 6.25, 8, 10, 2**32])
carect = np.array([0.579, 0.642, 0.69, 0.704, 0.757, 0.801, 0.830, 0.840, 0.872, 0.897, 0.917, 0.934, 0.947, 1])
rect_cafunc = inter(barect, carect, kind='quadratic')

#3d square prism coefficients as per RP H103
basq = np.array([1, 2, 3, 4, 5, 6, 7, 10])
casq = np.array([0.68, 0.36, 0.24, 0.19, 0.15, 0.13, 0.11, 0.08])
def sq_pris_cafunc(basq, a, b):
	return a*basq**b
	
popt, pcov = cf(sq_pris_cafunc, basq, casq)
	
#3d cylindrical coefficients as per RP H103
b2acyl = np.array([1.2, 2.5, 5, 9, 2**32])
cacyl = np.array([0.62, 0.78, 0.9, 0.96, 1])
try:
	cyl_cafunc = inter(b2acyl, cacyl, kind='cubic')
except np.linalg.linalg.LinAlgError:
	cyl_cafunc = inter(b2acyl, cacyl, kind='quadratic')	

'''
relevant DNV 4.6.3.2 The heave added mass of a surface piercing structure may be taken as half the heave added mass (in an
infinite fluid domain) of an object formed by the submerged part of the structure plus its mirror about the free surface.
'''

#2d ca function for rectangular bodies - surface piercing condition is included here only for future use. calling the addedmass_r2d() function will return an error string for a surface piercing body.
#However, this code could be used in future for surge and sway analysis if this program is extended, so it is being left in. Also note that a 2d circular ca function is not included - none of the
#platforms being modeled require this - if this is needed later for surge and sway analysis - adding this will be trivial.
def get_2d_rca(breadth, height, is_surface_piercing=False):
	if is_surface_piercing:
		return twod_cafunc(breadth/ (2*height)) / 2
	else:
		return twod_cafunc(breadth/height)

#Returns 3d add mass coefficients for rectangular plates
def get_rect_ca(b, a, is_surface_piercing=False):
	if is_surface_piercing:
		return (rect_cafunc(b/a))/2
	else:
		return rect_cafunc(b/a)
		
#3d ca coefficient for square prisms(in this case, the platform legs)	
def get_sq_pris_ca(height, breadth, is_surface_piercing=False):
	if is_surface_piercing:
		return (popt[0] * (2 * height / breadth)**popt[1])/2
	elif not is_surface_piercing:	
		return popt[0] * (height / breadth)**popt[1]	


#Returns heave added mass coefficient for circular disc/columns

def	get_circ_ca(is_surface_piercing=False):
	if is_surface_piercing:
		return 1/math.pi
	else:
		return 2/math.pi
	
def get_cyl_ca(height, radius, is_surface_piercing=False):
	try:
		if is_surface_piercing:
			return cyl_cafunc(height/radius)/2
		else:	
			return cyl_cafunc(height/(2*radius))
	except ValueError:
		if is_surface_piercing:
			return 1/math.pi
		else:
			return 2/math.pi
#now define reference volume functions for our basic shapes.

def get_2d_rar(breadth):
	return math.pi * (breadth/2)**2
	
def get_rect_vr(long_side, short_side):
	return (math.pi / 4)*(short_side**2)*long_side
	
def get_sq_pris_vr(height, breadth, is_surface_piercing):
	if is_surface_piercing:
		return breadth**2 * (2 * height)
	elif not is_surface_piercing:
		return breadth**2 * height
	
def get_circ_vr(radius):   
	return (4/3) * math.pi * radius**3
	
def get_cyl_vr(height, radius, is_surface_piercing=False):
	if get_cyl_ca(height, radius) == 2/math.pi:
		return (4/3) * math.pi * radius**3
	else:
		if is_surface_piercing:
			return math.pi * radius**2 *  (2 * height)
		elif not is_surface_piercing:
			return math.pi * radius**2 * height



#addedmass_r2d() function returns and added mass PER UNIT LENGTH - this must then be integrated or multiplied along the length of the body.

def addedmass_r2d(breadth, height, is_surface_piercing=False):
	if not is_surface_piercing:
		return rho * get_2d_rca(breadth, height, is_surface_piercing) * get_2d_rar(breadth)
	elif is_surface_piercing:
		print ('2d body added mass analysis is unsuitable for surface piercing bodies in heave - please use 3d analysis. surge and sway analysis will be added at a later date')  #should not be seen by users
		return 0		


def addedmass_rect_3d(b, a, height, is_surface_piercing=False):
	return rho * get_rect_ca(max(b, a),min(b, a), is_surface_piercing) * get_rect_vr(max(b, a), min(b, a))
	

def addedmass_sq_pris(height, breadth, is_surface_piercing=False):
	return rho * get_sq_pris_ca(height, breadth, is_surface_piercing) * get_sq_pris_vr(height, breadth, is_surface_piercing)
		
#these functions for circular and cylindrical body added mass were added after the decision to handle all surface piercing
#calculation in the subprocesses - thus these are not overdefined like the 2d and 3d rect functions above

def addedmass_circ_3d(radius, is_surface_piercing=False):
	return rho * get_circ_ca(is_surface_piercing) * get_circ_vr(radius)
	
def addedmass_cyl_3d(height, radius, is_surface_piercing=False):
	return rho * get_cyl_ca(height, radius, is_surface_piercing) * get_cyl_vr(height, radius)

'''
drag force components
'''	

plate_heave_drag_bh = [1, 5, 10, 2**32]
plate_heave_drag_cd = [1.16, 1.2, 1.5, 1.9]
plate_heave_cdfunc = inter(plate_heave_drag_bh, plate_heave_drag_cd, kind ='quadratic')

sq_pris_heave_drag_ld = [1, 1.5, 2, 2.5, 3, 4, 5]
sq_pris_heave_drag_cd = [1.15, 0.97, 0.87, 0.9, 0.93, 0.95, 0.95]
sq_pris_heave_drag_cdfunc = inter(sq_pris_heave_drag_ld, sq_pris_heave_drag_cd, kind = 'cubic')

cyl_heave_drag_ld = [0, 1, 2, 4, 7]
cyl_heave_drag_cd = [1.12, 0.91, 0.85, 0.87, 0.99]
cyl_heave_drag_cdfunc = inter(cyl_heave_drag_ld, cyl_heave_drag_cd, kind = 'cubic')

def get_plate_cd(length, breadth):
	return plate_heave_cdfunc(length/breadth)
	
def get_sq_pris_cd(length, square_side_length):
	return sq_pris_heave_drag_cdfunc(length/square_side_length)
	
def get_cyl_cd(length, diameter):
	return cyl_heave_drag_cdfunc(length/diameter)
	
'''
Wave heave motion functions, and dynamic pressure function
'''	
	
#Functions to define heave motions, i.e. z direction velocity and acceleration:
def z_vel():
	uz = ang_freq * w_amp * exp(w_num * z) * sin(w_num * x - ang_freq * t)
	return uz
	
def z_accel():
	az = -(ang_freq**2) * w_amp * exp(w_num * z) * cos(w_num * x - ang_freq * t)
	return az
	
def dyn_pres():
	p = rho * g * w_amp * exp(w_num*z)*cos(w_num * x - ang_freq * t)
	return p
	
'''
Inertia force functions
'''
	
# heave inertia force functions (suitable for heave only - need to update if surge and sway are added at a later date)

def heave_inertia_force_rect(breadth, height, length, xpos, zpos, is_surface_piercing=False):
	if not is_surface_piercing:
		try:
			ma = addedmass_r2d(breadth, height, is_surface_piercing)
			m2d = breadth * height * rho
			f = integrate((m2d+ma)*z_accel(), (x, xpos-length/2, xpos+length/2), conds='none')
		except ValueError:
			ma = addedmass_rect_3d(max(breadth, length), min(breadth, length), height, is_surface_piercing)
			f = ma * z_accel().subs(x, xpos) + integrate(breadth * dyn_pres(), (x, xpos-length/2, xpos+length/2), conds='none')
	elif is_surface_piercing:
		ma = addedmass_rect_3d(max(breadth, length), min(breadth, length), height, is_surface_piercing)
		f = ma * z_accel().subs(x, xpos) + integrate(breadth * dyn_pres(), (x, xpos-length/2, xpos+length/2), conds='none')
	return f.subs(z, zpos-(height/2))	

#assuming small diameter for now, this function is not suitabile for cylindrical platforms etc at the moment, (i.e. large diameters)
#as those would need to have angular integration. currently outside the scope of this project.
	
def heave_inertia_force_circ(radius, height, xpos, zpos, is_surface_piercing=False):
	ma = addedmass_circ_3d(radius, is_surface_piercing)
	m = math.pi * radius**2 * height * rho
	return (m + ma) * z_accel().subs(x, xpos).subs(z, zpos)
	
def heave_inertia_force_cyl(height, radius, xpos, zpos, is_surface_piercing=False):
	ma = addedmass_cyl_3d(height, radius, is_surface_piercing)
	m = math.pi * radius**2 * height * rho
	return (m + ma) * z_accel().subs(x, xpos).subs(z, zpos)
	
'''
Drag Force functions
'''

def plate_drag_force(length, breadth, xpos, zpos):
	f = (1/2) * rho * get_plate_cd(max(length, breadth), min(length, breadth)) * breadth * z_vel()**2
	posneg_check = integrate(z_vel(), (x, xpos-length/2, xpos+length/2), conds='none')
	f_n = integrate(f, (x, xpos-length/2, xpos+length/2), conds='none')
	if posneg_check.subs(z, zpos).subs(t, 0).subs(w_num, 1).subs(ang_freq, 1) > 0:
		return expand_trig(simplify(f_n.subs(z, zpos))).subs(t, 0)
	else:
		return expand_trig(simplify(-f_n.subs(z, zpos))).subs(t, 0)
		
def sq_pris_drag_force(length, square_side_length, xpos, zpos):
	f = (1/2) * rho * get_sq_pris_cd(length, square_side_length) * square_side_length * z_vel()**2
	posneg_check = integrate(z_vel(), (x, xpos-square_side_length/2, xpos+square_side_length/2), conds='none')
	f_n = integrate(f, (x, xpos-square_side_length/2, xpos+square_side_length/2), conds='none')
	if posneg_check.subs(z, zpos).subs(t, 0).subs(w_num, 1).subs(ang_freq, 1) > 0:
		return expand_trig(simplify(f_n.subs(z, zpos-length/2))).subs(t, 0)
	else:
		return expand_trig(simplify(-f_n.subs(z, zpos-length/2))).subs(t, 0)

#this assumes small diameter - unsuitable for cylindrical platforms etc - may be added in future.		
		
def cyl_drag_force(length, diameter, xpos, zpos):
	f = (1/2) * rho * get_cyl_cd(length, diameter) * (math.pi/4) * diameter**2 * z_vel()**2
	if z_vel().subs(z, zpos).subs(x, xpos).subs(t, 0).subs(w_num, 1).subs(ang_freq, 1) > 0:
		return expand_trig(simplify(f.subs(z, zpos-length/2))).subs(t, 0)
	else:
		return expand_trig(simplify(-f.subs(z, zpos-length/2))).subs(t, 0)		
	
	
'''
Platform Model Functions
'''

def simp_platform_force(draft, square_leg_width, length, pontoon_length, pontoon_width, pontoon_height):
	simp_platform_DF = ((sq_pris_drag_force(draft, square_leg_width, -((length/2)-(square_leg_width/2)), -draft/2) 			#square leg(s) in negative x direction
		+ sq_pris_drag_force(draft, square_leg_width, ((length/2)-(square_leg_width/2)), -draft/2)) * 2 	  				#square legs in pos x direction
		+ 2 * (plate_drag_force(pontoon_length, pontoon_width, 0, -(draft-(pontoon_height/2)))) 				  			#2 pontoons in x direction lengthwise
		+ plate_drag_force(pontoon_width, pontoon_length, (square_leg_width/2)-length/2, -(draft-(pontoon_height/2)))  		# rear pontoon
		+ plate_drag_force(pontoon_width, pontoon_length, -(square_leg_width/2)+length/2, -(draft-(pontoon_height/2)))) 	#front pontoon	
	simp_platform_IF = ((heave_inertia_force_rect(square_leg_width, draft, square_leg_width, -((length/2)-(square_leg_width/2)), -draft/2,  True) 	#square leg(s) in negative x direction
		+ heave_inertia_force_rect(square_leg_width, draft, square_leg_width, ((length/2)-(square_leg_width/2)), -draft/2, True)) * 2 	  			#square legs in pos x direction
		+ 2 * (heave_inertia_force_rect(pontoon_width, pontoon_height, pontoon_length, 0, -(draft-(pontoon_height/2)))) 				  			#2 pontoons in x direction lengthwise
		+ heave_inertia_force_rect(pontoon_length, pontoon_height, pontoon_width, (square_leg_width/2)-length/2, -(draft-(pontoon_height/2)))  		# rear pontoon
		+ heave_inertia_force_rect(pontoon_length, pontoon_height, pontoon_width, -(square_leg_width/2)+length/2, -(draft-(pontoon_height/2)))) 	#front pontoon
	return simp_platform_DF + simp_platform_IF
	
def plat_with_heave_plate_force(heave_plate_separation, draft, length, heave_plate_cutout_width, heave_plate_height, heave_plate_section_width, extensor_diameter, extensor_sep, extensor_spacing):		
	plat_wheave_pl_DF = ( 
		 2 * plate_drag_force(length, heave_plate_section_width, 0, -(draft + heave_plate_separation + heave_plate_height/2))																# 2 longitudinal sections
		+ plate_drag_force(heave_plate_section_width, heave_plate_cutout_width, -(length - heave_plate_section_width)/2, -(draft + heave_plate_separation + heave_plate_height/2))			# aft small section
		+ plate_drag_force(heave_plate_section_width, heave_plate_cutout_width, (length - heave_plate_section_width)/2, -(draft + heave_plate_separation + heave_plate_height/2)))			#fore small zection
	plat_wheave_pl_IF = (
		2 * (heave_inertia_force_rect(heave_plate_section_width, heave_plate_height, length-(2*heave_plate_section_width), 0, -(draft+heave_plate_separation+heave_plate_height/2))							#2 pontoons in x direction lengthwise
		+ heave_inertia_force_rect(heave_plate_section_width, heave_plate_height, heave_plate_section_width, -(length-heave_plate_section_width)/2, -(draft+heave_plate_separation+heave_plate_height/2))	#aft square sections
		+ heave_inertia_force_rect(heave_plate_section_width, heave_plate_height, heave_plate_section_width, (length-heave_plate_section_width)/2, -(draft+heave_plate_separation+heave_plate_height/2)))	#fore square sections
		+ heave_inertia_force_rect(heave_plate_cutout_width, heave_plate_height, heave_plate_section_width, (heave_plate_section_width/2)-length/2, -(draft+heave_plate_separation+heave_plate_height/2))  	# aft heave plate section
		+ heave_inertia_force_rect(heave_plate_cutout_width, heave_plate_height, heave_plate_section_width, -(heave_plate_section_width/2)+length/2, -(draft+heave_plate_separation+heave_plate_height/2)) 	# fore heave plate section	
		- 2 * (heave_inertia_force_cyl(heave_plate_height, extensor_diameter/2, -extensor_spacing/2, -(draft + heave_plate_separation + heave_plate_height/2)) 								#outer rear extensors
		+ heave_inertia_force_cyl(heave_plate_height, extensor_diameter/2, -((extensor_spacing/2)+(extensor_sep/math.sqrt(2))), -(draft + heave_plate_separation + heave_plate_height/2))	#inner rear extensors
		+ heave_inertia_force_cyl(heave_plate_height, extensor_diameter/2, extensor_spacing/2, -(draft + heave_plate_separation + heave_plate_height/2))									#outer fore extensors 
		+ heave_inertia_force_cyl(heave_plate_height, extensor_diameter/2, (extensor_spacing/2)+(extensor_sep/math.sqrt(2)), -(draft + heave_plate_separation + heave_plate_height/2))))	#inner fore extensors
	return plat_wheave_pl_DF + plat_wheave_pl_IF + simp_platform_force(eds_draft, eds_square_leg_width, eds_length, eds_pontoon_length, eds_pontoon_width, eds_pontoon_height)	

'''
Function to perform DDS Analysis
'''

if dds_do:	
	simp_platform_FTOTAL = simp_platform_force(dds_draft, dds_square_leg_width, dds_length, dds_pontoon_length, dds_pontoon_width, dds_pontoon_height)
	simp_platform_AM = 	((4 * addedmass_rect_3d(dds_square_leg_width, dds_square_leg_width, dds_draft,  True) 																# 4 square legs
		+ 2 * (addedmass_r2d(dds_pontoon_width, dds_pontoon_height)*dds_pontoon_length + addedmass_r2d(dds_pontoon_length, dds_pontoon_height) * dds_pontoon_width)))		#4 pontoons
	simp_platform_MASS = 4 * rho * (dds_square_leg_width**2 * dds_draft + dds_pontoon_height * dds_pontoon_length * dds_pontoon_width)
	simp_plat_k33 = rho * g * dds_square_leg_width**2 * 4
	simp_plat_wn = math.sqrt(simp_plat_k33/(simp_platform_AM + simp_platform_MASS))
	
	ifv_s = [simp_platform_FTOTAL.subs(w_num, x).subs(ang_freq, y).subs(t, 0) for x, y in zip(k, w_l)]
	cap_omega_s = [ x / simp_plat_wn for x in w_l]
	denom_s = [ simp_plat_k33 * math.sqrt((1-x**2)**2 +(2 * x * zeta)**2) for x in cap_omega_s]
	x3_s = [abs(x/y) for x, y in zip(ifv_s, denom_s)]
	response_spec_s = [x * y**2 for x, y in zip(spec, x3_s)]
	print('Processing Deep Draft Semi')
	if ret_mass:
		print ('DDS Mass = {:.2f} Kg (assuming density in Kg/m^3)'.format(simp_platform_MASS))
	if ret_addedmass:
		print ('DDS Added Mass = {:.2f} Kg (assuming density in Kg/m^3)'.format(simp_platform_AM))
	if ret_wn:
		print ('DDS Natural Frequency = {:.4f} rad/sec'.format(simp_plat_wn))
	if ret_tn:
		print ('DDS Natural Period = {:.2f} sec'.format(2*math.pi/simp_plat_wn))
	print('+----------------------------------------+')	
	if ret_sig or ret_mean or ret_max:
		print('Calculating Characteristic Values')
		m_r_s = scipy.integrate.simps(response_spec_s, w_l)
		if ret_sig:
			print('DDS Significant Response Amplitude = {:.4f}m'.format(2 * math.sqrt(m_r_s)))
		if ret_mean:
			print('DDS Mean Response Amplitude = {:.4f}m'.format(1.25 * math.sqrt(m_r_s)))
		if ret_max:
			print('DDS Maximum Response Amplitude = {:.4f}m'.format(4 * math.sqrt(m_r_s)))
		print('+----------------------------------------+')	
		
'''
Function to perform EDS Analysis
'''
		
if eds_do:	
	plat_wheave_pl_FTOTAL = plat_with_heave_plate_force(heave_plate_separation, eds_draft, eds_length, heave_plate_cutout_width, heave_plate_height, heave_plate_section_width, extensor_diameter, extensor_sep, extensor_spacing)	
	plat_wheave_pl_AM = ((4 * addedmass_rect_3d(eds_square_leg_width, eds_square_leg_width, eds_draft,  True) 
		+ 2 * (addedmass_r2d(eds_pontoon_width, eds_pontoon_height)*eds_pontoon_length + addedmass_r2d(eds_pontoon_length, eds_pontoon_height) * eds_pontoon_width)) +
		4 * addedmass_r2d(eds_square_leg_width, heave_plate_height) * (eds_length-heave_plate_section_width)
		- 8 * addedmass_circ_3d(extensor_diameter/2, True))	
	plat_wheave_pl_MASS = ((4 * rho * (eds_square_leg_width**2 * eds_draft + eds_pontoon_height * eds_pontoon_length * eds_pontoon_width))
		+ rho * (heave_plate_height * ( eds_length**2 - heave_plate_cutout_width**2 )+ 8 * math.pi * (extensor_diameter**2 / 4) * heave_plate_separation))	
	plat_wheave_pl_k33 = rho * g * eds_square_leg_width**2 * 4
	plat_wheave_pl_wn = math.sqrt(plat_wheave_pl_k33/(plat_wheave_pl_AM + plat_wheave_pl_MASS))	
	
	ifv = [plat_wheave_pl_FTOTAL.subs(w_num, x).subs(ang_freq, y).subs(t, 0) for x, y in zip(k, w_l)]
	cap_omega = [ x / plat_wheave_pl_wn for x in w_l]
	denom = [ plat_wheave_pl_k33 * math.sqrt((1-x**2)**2 +(2 * x * zeta)**2) for x in cap_omega]
	x3 = [abs(x/y) for x, y in zip(ifv, denom)]
	response_spec = [x * y**2 for x, y in zip(spec, x3)]
	print('Processing Extendable Draft Semi')
	if ret_mass:
		print ('EDS Mass = {:.2f} Kg (assuming density in Kg/m^3)'.format(plat_wheave_pl_MASS))
	if ret_addedmass:
		print ('EDS Added Mass = {:.2f} Kg (assuming density in Kg/m^3)'.format(plat_wheave_pl_AM))		
	if ret_wn:
		print ('EDS Natural Frequency = {:.4f} rad/sec'.format(plat_wheave_pl_wn))
	if ret_tn:
		print ('EDS Natural Period = {:.2f} sec'.format(2*math.pi/plat_wheave_pl_wn))
	print('+----------------------------------------+')		
	if ret_sig or ret_mean or ret_max:
		print('Calculating Characteristic Values')
		m_r = scipy.integrate.simps(response_spec, w_l)
		if ret_sig:
			print('EDS Significant Response Amplitude = {:.4f}m'.format(2 * math.sqrt(m_r)))
		if ret_mean:
			print('EDS Mean Response Amplitude = {:.4f}m'.format(1.25 * math.sqrt(m_r)))
		if ret_max:
			print('EDS Maximum Response Amplitude = {:.4f}m'.format(4 * math.sqrt(m_r)))
		print('+----------------------------------------+')		

'''
Function to perform Varying Depth Analysis
'''		

if var_do:
	print('Processing Geometry for Varying Heave Plate Depth')
	if not eds_do:
		plat_wheave_pl_AM = ((4 * addedmass_rect_3d(eds_square_leg_width, eds_square_leg_width, eds_draft,  True) 
		+ 2 * (addedmass_r2d(eds_pontoon_width, eds_pontoon_height)*eds_pontoon_length + addedmass_r2d(eds_pontoon_length, eds_pontoon_height) * eds_pontoon_width)) +
		4 * addedmass_r2d(eds_square_leg_width, heave_plate_height) * (eds_length-heave_plate_section_width)
		- 8 * addedmass_circ_3d(extensor_diameter/2, True))	
		plat_wheave_pl_MASS = ((4 * rho * (eds_square_leg_width**2 * eds_draft + eds_pontoon_height * eds_pontoon_length * eds_pontoon_width))
		+ rho * (heave_plate_height * ( eds_length**2 - heave_plate_cutout_width**2 )+ 8 * math.pi * (extensor_diameter**2 / 4) * heave_plate_separation))	
		plat_wheave_pl_k33 = rho * g * eds_square_leg_width**2 * 4
		plat_wheave_pl_wn = math.sqrt(plat_wheave_pl_k33/(plat_wheave_pl_AM + plat_wheave_pl_MASS))	

	hp_sep_vals = [x for x in frange(sep_s, sep_e, sep_divs+1)]
	cap_omega = [ x / plat_wheave_pl_wn for x in w_l]
	denom = [ plat_wheave_pl_k33 * math.sqrt((1-x**2)**2 +(2 * x * zeta)**2) for x in cap_omega]
	ifv_var = []
	x3_var = []
	response_spec_var = []
	print('Done')
	if ret_RAO or ret_response:
		print('Processing Results for Each Heave Plate Depth')
		for index, value in enumerate(hp_sep_vals):
			plat_wheave_pl_FTOTAL = plat_with_heave_plate_force(value, eds_draft, eds_length, heave_plate_cutout_width, heave_plate_height, heave_plate_section_width, extensor_diameter, extensor_sep, extensor_spacing)
			ifv_var.append([plat_wheave_pl_FTOTAL.subs(w_num, x).subs(ang_freq, y).subs(t, 0) for x, y in zip(k, w_l)])			
			if ret_RAO:
				x3_var.append([abs(x/y) for x, y in zip(ifv_var[index], denom)])
			if ret_response:
				response_spec_var.append([x * y**2 for x, y in zip(spec, x3_var[index])])
				if ret_sig or ret_mean or ret_max:
					m_r_var = scipy.integrate.simps(response_spec_var[index], w_l)
					if ret_sig:
						print('EDS Significant Response Amplitude, {:.2f}m Plate Seperation = {:.4f}m'.format(hp_sep_vals[index], 2 * math.sqrt(m_r_var)))
					if ret_mean:
						print('EDS Mean Response Amplitude, {:.2f}m Plate Seperation = {:.4f}m'.format(hp_sep_vals[index], 1.25 * math.sqrt(m_r_var)))
					if ret_max:
						print('EDS Maximum Response Amplitude, {:.2f}m Plate Seperation = {:.4f}m'.format(hp_sep_vals[index], 4 * math.sqrt(m_r_var)))	
					print('+----------------------------------------+')		
	else:
		print('Varying heave plate depth analysis has been selected, but no results requested. Please note, both RAO and Response Spectrum are needed in order to calculate characteristic response values. Please enable these options to complete the analysis')			

'''
Results Export to Excel
'''
		
if ret_xls:
	print('Exporting Results to Excel...')
	sfw = Series(w_l, name='Frequency (rad/sec)')
	sp = Series(period, name='Period (s)')
	series_list = [sfw,sp]
	if ret_spec:
		if pms_do:
			s1 = Series([round(x,10) for x in spec], name='Sea State {} Spectrum'.format(int(sea_state)))
			series_list.append(s1)	
		elif not pms_do:
			s1 = Series([round(x,10) for x in spec], name='User Defined Spectrum')
			series_list.append(s1)	
	if dds_do:
		if ret_RAO:
			s2 = Series([round(x,10) for x in x3_s], name='DDS Heave RAO (m/m)')
			series_list.append(s2)	
		if ret_response:
			s3 = Series([round(x,10) for x in response_spec_s], name='DDS Response')
			series_list.append(s3)		
	if eds_do:
		if ret_RAO:
			s4 = Series([round(x,10) for x in x3], name='EDS Heave RAO (m/m)')
			series_list.append(s4)	
		if ret_response:
			s5 = Series([round(x,10) for x in response_spec], name='EDS Response')
			series_list.append(s5)
	if var_do:
		if ret_RAO:
			for index, value in enumerate(x3_var):
				series_list.append(Series([round(x,10) for x in value], name='EDS RAO, {:.2f}m Plate Seperation'.format(hp_sep_vals[index])))
		if ret_response:
			for index, value in enumerate(response_spec_var):
				series_list.append(Series([round(x,10) for x in value], name='EDS Response, {:.2f}m Plate Seperation'.format(hp_sep_vals[index])))
	df = pandas.concat(series_list, axis=1)
	df.to_excel('Heavecalc_Results_{}.xlsx'.format(dt.datetime.strftime(dt.datetime.now(), '%d%m%Y_%H.%M.%S')), sheet_name='Heavecalc_Results', index=False)
	print('Export Complete - Excel results file is located in the working directory')	

'''
Generate result Plots
'''	
		
if ret_plot:
	print('Plotting Selected Results')
	if ret_x_w:
		if ret_spec:
			if pms_do:
				plot.plot(w_l, spec, 'b:', label='Sea State {} Spectrum'.format(int(sea_state)))
			elif not pms_do:
				plot.plot(w_l, spec, 'b:', label='User Defined Spectrum')
			
		if dds_do:
			if ret_RAO:
				plot.plot(w_l, x3_s, 'b-', label='DDS Heave RAO')
			if ret_response:
				plot.plot(w_l, response_spec_s, 'y--', label='DDS Response')
		
		if eds_do:
			if ret_RAO:
				plot.plot(w_l, x3, 'g-', label='EDS Heave RAO')
			if ret_response:
				plot.plot(w_l, response_spec, 'r--', label='EDS Response')
		
		if var_do:
			if ret_RAO:
				for index, value in enumerate(x3_var):
					plot.plot(w_l, value, label='EDS RAO, {:.2f}m Plate Seperation'.format(hp_sep_vals[index]))
			if ret_response:		
				for index, value in enumerate(response_spec_var):
					plot.plot(w_l, value, label='EDS Response, {:.2f}m Plate Seperation'.format(hp_sep_vals[index]))
		plot.xlabel('Frequency (rad/s)')	
		
	elif not ret_x_w:
		if ret_spec:
			if pms_do:
				plot.plot(period, spec, 'b:', label='Sea State {} Spectrum'.format(int(sea_state)))
			elif not pms_do:
				plot.plot(period, spec, 'b:', label='User Defined Spectrum')
			
		if dds_do:
			if ret_RAO:
				plot.plot(period, x3_s, 'b-', label='DDS Heave RAO')
			if ret_response:
				plot.plot(period, response_spec_s, 'y--', label='DDS Response')
		
		if eds_do:
			if ret_RAO:
				plot.plot(period, x3, 'g-', label='EDS Heave RAO')
			if ret_response:
				plot.plot(period, response_spec, 'r--', label='EDS Response')
				
		if var_do:
			if ret_RAO:
				for index, value in enumerate(x3_var):
					plot.plot(period, value, label='EDS RAO, {:.2f}m Plate Seperation'.format(hp_sep_vals[index]))
			if ret_response:		
				for index, value in enumerate(response_spec_var):
					plot.plot(period, value, label='EDS Response, {:.2f}m Plate Seperation'.format(hp_sep_vals[index]))				
		plot.xlabel('Wave Period (s)')	
	
	plot.ylabel('RAO(m/m)')	
	plot.legend(loc=1)
	plot.show()	
	
'''
old test code below this point,
not accessible to users through the settings file. May be integrated to the main code in future if required.
can be accessed if needed by setting all return options in the configuration to 0, and manually setting the required
code below to True.
'''			

#Plot of force on EDS platform across frequency range.
if False:
	plat_wheave_pl_FTOTAL = plat_with_heave_plate_force(heave_plate_separation, eds_draft, eds_length, heave_plate_cutout_width, heave_plate_height, heave_plate_section_width, extensor_diameter, extensor_sep, extensor_spacing)
	ifv = [plat_wheave_pl_FTOTAL.subs(w_num, x).subs(ang_freq, y).subs(t, 0) for x, y in zip(k, w_l)]
	plot.plot(w_l, ifv, 'g-')
	plot.show()	

input('Program complete, please hit [return] key to exit.')
