text: from openmdao.api import ExplicitComponent
  import numpy as np
  import os
  import sys
  from wisdem.pymap import pyMAP
  
  from wisdem.commonse import gravity, Enum
  from wisdem.commonse.utilities import assembleI, unassembleI
  
  Anchor    = Enum('DRAGEMBEDMENT SUCTIONPILE')
  NLINES_MAX = 15
  NPTS_PLOT = 20
  
  class MapMooring(ExplicitComponent):
      """
      OpenMDAO Component class for mooring system attached to sub-structure of floating offshore wind turbines.
      Should be tightly coupled with Spar class for full system representation.
      """
  
      def setup(self):
      
          # Variables local to the class and not OpenMDAO
          self.min_break_load      = None
          self.wet_mass_per_length = None
          self.axial_stiffness     = None
          self.area                = None
          self.cost_per_length     = None
          self.finput              = None
          self.tlpFlag             = False
  
          # Environment
          self.add_input('water_density', val=0.0, units='kg/m**3', desc='density of water')
          self.add_input('water_depth', val=0.0, units='m', desc='water depth')
  
          # Material properties
  
          # Inputs from SparGeometry
          self.add_input('fairlead_radius', val=0.0, units='m', desc='Outer spar radius at fairlead depth (point of mooring attachment)')
          
          # Design variables
          self.add_input('fairlead', val=0.0, units='m', desc='Depth below water for mooring line attachment')
          self.add_input('mooring_line_length', val=0.0, units='m',desc='Unstretched total mooring line length')
          self.add_input('anchor_radius', val=0.0, units='m', desc='radius from center of spar to mooring anchor point')
          self.add_input('mooring_diameter', val=0.0, units='m',desc='diameter of mooring line')
  
          # User inputs (could be design variables)
          self.add_input('number_of_mooring_connections', val=3, desc='number of mooring connections on vessel')
          self.add_input('mooring_lines_per_connection', val=1, desc='number of mooring lines per connection')
          self.add_discrete_input('mooring_type', val='CHAIN', desc='chain, nylon, polyester, fiber, or iwrc')
          self.add_discrete_input('anchor_type', val='DRAGEMBEDMENT', desc='SUCTIONPILE or DRAGEMBEDMENT')
          self.add_input('max_offset', val=0.0, units='m',desc='X offsets in discretization')
          self.add_input('operational_heel', val=0.0, units='deg',desc='Maximum angle of heel allowable during operation')
          self.add_input('max_survival_heel', val=0.0, units='deg', desc='max heel angle for turbine survival')
          self.add_input('gamma_f', val=0.0, desc='Safety factor for mooring line tension')
  
          # Cost rates
          self.add_input('mooring_cost_factor', val=0.0, desc='miscellaneous cost factor in percent')
  
          # Outputs
          self.add_output('number_of_mooring_lines', val=0, desc='total number of mooring lines')
          self.add_output('mooring_mass', val=0.0, units='kg',desc='total mass of mooring')
          self.add_output('mooring_moments_of_inertia', val=np.zeros(6), units='kg*m**2', desc='mass moment of inertia of mooring system about fairlead-centerline point [xx yy zz xy xz yz]')
          self.add_output('mooring_cost', val=0.0, units='USD',desc='total cost for anchor + legs + miscellaneous costs')
          self.add_output('mooring_stiffness', val=np.zeros((6,6)), units='N/m', desc='Linearized stiffness matrix of mooring system at neutral (no offset) conditions.')
          self.add_output('anchor_cost', val=0.0, units='USD',desc='total cost for anchor')
          self.add_output('mooring_neutral_load', val=np.zeros((NLINES_MAX,3)), units='N',desc='mooring vertical load in all mooring lines')
          self.add_output('max_offset_restoring_force', val=0.0, units='N',desc='sum of forces in x direction after max offset')
          self.add_output('operational_heel_restoring_force', val=np.zeros((NLINES_MAX,3)), units='N',desc='forces for all mooring lines after operational heel')
          self.add_output('mooring_plot_matrix', val=np.zeros((NLINES_MAX, NPTS_PLOT, 3)), units='m', desc='data matrix for plotting') 
  
          # Output constriants
          self.add_output('axial_unity', val=0.0, units='m',desc='range of damaged mooring')
          self.add_output('mooring_length_max', val=0.0, desc='mooring line length ratio to nodal distance')
          
          # Derivatives
          self.declare_partials('*', '*', method='fd', form='central', step=1e-6)
  
          
      def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
          """Sets mooring line properties then writes MAP input file and executes MAP.
          
          INPUTS:
          ----------
          inputs   : dictionary of input parameters
          outputs : dictionary of output parameters
          
          OUTPUTS  : none (all unknown dictionary values set)
          """
  
          # Set characteristics based on regressions / empirical data
          self.set_properties(inputs, discrete_inputs)
  
          # Set geometry profile
          self.set_geometry(inputs, outputs)
  
          # Write MAP input file and analyze the system at every angle
          self.runMAP(inputs, discrete_inputs, outputs)
  
          # Compute costs for the system
          self.compute_cost(inputs, discrete_inputs, outputs)
  
  
      def set_properties(self, inputs, discrete_inputs):
          """Sets mooring line properties: Minimum Breaking Load, Mass per Length, Axial Stiffness, Cross-Sectional Area, Cost-per-Length.
          
          INPUTS:
          ----------
          inputs   : dictionary of input parameters
          
          OUTPUTS  : Parameters are class variables and are set internally
  
          References:
          https://daim.idi.ntnu.no/masteroppgaver/015/15116/masteroppgave.pdf
          http://offshoremechanics.asmedigitalcollection.asme.org/article.aspx?articleid=2543338
          https://www.orcina.com/SoftwareProducts/OrcaFlex/Documentation/Help/Content/html/
          Chain.htm
          Chain,AxialandBendingStiffness.htm
          Chain,MechanicalProperties.htm
          RopeWire.htm
          RopeWire,MinimumBreakingLoads.htm
          RopeWire,Massperunitlength.htm
          RopeWire,AxialandBendingStiffness.htm
          """
          
          # Unpack variables
          Dmooring       = inputs['mooring_diameter']
          lineType       = discrete_inputs['mooring_type'].upper()
  
          # Set parameters based on regressions for different mooring line type
          Dmooring2 = Dmooring**2
  
          # TODO: Costs per unit length are not synced with new input sources
          if lineType == 'CHAIN':
              self.min_break_load      = 2.74e7  * Dmooring2 * (44.0 - 80.0*Dmooring)
              # Use a linear fit to the other fit becuase it is poorly conditioned for optimization
              #self.min_break_load      = 1e3*np.maximum(1.0, -5445.2957034820683+176972.68498888266*Dmooring)
              self.wet_mass_per_length = 19.9e3  * Dmooring2 # From Orca, 7983.34117 OC3 definiton doc
              self.axial_stiffness     = 8.54e10 * Dmooring2 # From Orca, 4.74374e10 OC3 definiton doc, 
              self.area                = 2.0 * 0.25 * np.pi * Dmooring2
              self.cost_per_length     = 3.415e4  * Dmooring2 #0.58*1e-3*self.min_break_load/gravity - 87.6
  
          elif lineType == 'NYLON':
              self.min_break_load      = 139357e3 * Dmooring2
              self.wet_mass_per_length = 0.6476e3 * Dmooring2
              self.axial_stiffness     = 1.18e8   * Dmooring2
              self.area                = 0.25 * np.pi * Dmooring2
              self.cost_per_length     = 3.415e4  * Dmooring2 #0.42059603*1e-3*self.min_break_load/gravity + 109.5
  
          elif lineType == 'POLYESTER':
              self.min_break_load      = 170466e3 * Dmooring2
              self.wet_mass_per_length = 0.7978e3 * Dmooring2
              self.axial_stiffness     = 1.09e9   * Dmooring2
              self.area                = 0.25 * np.pi * Dmooring2
              self.cost_per_length     = 3.415e4  * Dmooring2 #0.42059603*1e-3*self.min_break_load/gravity + 109.5
  
          elif lineType == 'FIBER': # Wire rope with fiber rope
              self.min_break_load      = 584175e3 * Dmooring2
              self.wet_mass_per_length = 3.6109e3 * Dmooring2
              self.axial_stiffness     = 3.67e10  * Dmooring2
              self.area                = 0.455 * 0.25 * np.pi * Dmooring2
              self.cost_per_length     = 2.0 * 6.32e4  * Dmooring2 #0.53676471*1e-3*self.min_break_load/gravity
  
          elif lineType == 'IWRC': # Wire rope with steel core
              self.min_break_load      = 633358e3 * Dmooring2
              self.wet_mass_per_length = 3.9897e3 * Dmooring2
              self.axial_stiffness     = 4.04e10  * Dmooring2
              self.area                = 0.455 * 0.25 * np.pi * Dmooring2
              self.cost_per_length     = 6.32e4  * Dmooring2 #0.33*1e-3*self.min_break_load/gravity + 139.5
  
          else:
              raise ValueError('Available line types are: chain nylon polyester fiber iwrc')
  
              
      def set_geometry(self, inputs, outputs):
          # Unpack variables
          fairleadDepth = inputs['fairlead']
          R_fairlead    = inputs['fairlead_radius']
          R_anchor      = inputs['anchor_radius']
          waterDepth    = inputs['water_depth']
          L_mooring     = inputs['mooring_line_length']
          max_heel      = inputs['max_survival_heel']
          gamma         = inputs['gamma_f']
          
          if L_mooring > (waterDepth - fairleadDepth):
              self.tlpFlag = False
              
              # Create constraint that line isn't too long that there is no catenary hang
              outputs['mooring_length_max'] = L_mooring / (0.95 * (R_anchor + waterDepth - fairleadDepth) )
          else:
              self.tlpFlag = True
              # Create constraint that we don't lose line tension
              outputs['mooring_length_max'] = L_mooring / ( (waterDepth - fairleadDepth - gamma*R_fairlead*np.sin(np.deg2rad(max_heel))) )
  
      
      def write_line_dictionary(self, inputs, discrete_inputs, cable_sea_friction_coefficient=0.65):
          """Writes LINE DICTIONARY section of input.map file
          
          INPUTS:
          ----------
          inputs   : dictionary of input parameters
          cable_sea_friction_coefficient : coefficient of friction with sea floor (defaults to 0.65)
          
          OUTPUTS  : none
          """
          # Unpack variables
          rhoWater = inputs['water_density']
          lineType = discrete_inputs['mooring_type'].lower()
          Dmooring = inputs['mooring_diameter']
          
          self.finput.append('---------------------- LINE DICTIONARY ---------------------------------------')
          self.finput.append('LineType  Diam      MassDenInAir   EA            CB   CIntDamp  Ca   Cdn    Cdt')
          self.finput.append('(-)       (m)       (kg/m)        (N)           (-)   (Pa-s)    (-)  (-)    (-)')
          self.finput.append('%s   %.5f   %.5f   %.5f   %.5f   1.0E8   0.6   -1.0   0.05' %
                            (lineType, Dmooring, self.wet_mass_per_length, self.axial_stiffness, cable_sea_friction_coefficient) )
  
          
      def write_node_properties_header(self):
          """Writes NODE PROPERTIES section header of input.map file
          
          INPUTS: none 
          ----------
          
          OUTPUTS  : none
          """
          self.finput.append('---------------------- NODE PROPERTIES ---------------------------------------')
          # Doesn't like some weird character here somewhere
          #self.finput.append('Node  Type       X       Y       Z      M     B     FX      FY      FZ')
          #self.finput.append('(-)   (-)       (m)     (m)     (m)    (kg)  (m^3)  (N)     (N)     (N)')
          self.finput.append('Node Type X     Y    Z   M     V FX FY FZ')
          self.finput.append('(-)  (-) (m)   (m)  (m) (kg) (m^3) (kN) (kN) (kN)')
  
  
  
      def write_node_properties(self, number, node_type, x_pos, y_pos, z_pos,
                                point_mass=0, displaced_volume=0,
                                x_force=None, y_force=None, z_force=None):
          """Writes NODE PROPERTIES data of input.map file.  Nodes are connections between mooring lines and bridles, vessels, and anchors
          
          INPUTS:
          ----------
          number           : The node number listing in the input file
          node_type        : fix / vessel / connect (fix=anchor)
          x_, y_, z_pos    : position of node in coordinate system (separate inputs)
          point_mass       : see MAP reference (defaults to 0)
          displaced_volume : see MAP reference (defaults to 0)
          x_, y_, z_force  : see MAP reference (defaults to None)
          
          OUTPUTS  : none
          """
          # Ensure this connection is something MAP understands
          nodeStr = node_type.lower()
          if not nodeStr in ['fix', 'connect', 'vessel']:
              raise ValueError('%s is not a valid node type for node %d' % (node_type, number))
  
          # If this is a node between two lines have to specify connection details
          if nodeStr == 'connect':
              try:
                  x_force = float(x_force)
                  y_force = float(y_force)
                  z_force = float(z_force)
              except:
                  raise ValueError('%s must have numerical force applied values.' % node_type)
  
          # Set location strings
          forceStr = '#   #   #' 
          if nodeStr == 'connect':
              forceStr = '%f   %f   %f' % (x_force, y_force, z_force)
              posStr   = '#%f   #%f   #%f   ' % (x_pos, y_pos, z_pos)
          elif nodeStr == 'fix':
              posStr   = '%f   %f   depth   ' % (x_pos, y_pos)
          elif nodeStr == 'vessel':
              posStr   = '%f   %f   %f   ' % (x_pos, y_pos, z_pos)
  
          # Write the connection line
          line = ('%d   ' % number)
          line += ('%s   ' % node_type)
          line += (posStr)
          line += ('%f   %f   ' % (point_mass, displaced_volume) )
          line += (forceStr)
          self.finput.append(line)
  
      def write_line_properties(self, inputs, discrete_inputs, line_number=1, anchor_node=2, fairlead_node=1, flags=''):
          """Writes LINE PROPERTIES section of input.map file that connects multiple nodes
          
          INPUTS:
          ----------
          inputs        : dictionary of input parameters- only 'mooring_type' is used
          line_number   : Line ID number (defaults to 1)
          anchor_node   : Node number corresponding to anchor (defaults to 1)
          fairlead_node : Node number corresponding to fairlead (vessel) node (defaults to 2)
          flags         : see MAP reference (defaults to empty string ' ')
          
          OUTPUTS  : none
          """
          # Add flag for taut lines
          if self.tlpFlag:
              flags += ' LINEAR SPRING'
              
          self.finput.append('---------------------- LINE PROPERTIES ---------------------------------------')
          self.finput.append('Line    LineType  UnstrLen  NodeAnch  NodeFair  Flags')
          self.finput.append('(-)      (-)       (m)       (-)       (-)       (-)')
          self.finput.append('%d   %s   %f   %d   %d   %s' %
                            (line_number, discrete_inputs['mooring_type'], inputs['mooring_line_length'], anchor_node, fairlead_node, flags) )
  
          
      def write_solver_options(self, inputs):
          """Writes SOLVER OPTIONS section of input.map file, 
          which includes repeating node/line arrangement in even angular spacing around structure.
          
          INPUTS:
          ----------
          inputs        : dictionary of input parameters
          
          OUTPUTS  : none
          """
          # Unpack variables
          n_connect     = max(1, int(inputs['number_of_mooring_connections']))
          
          self.finput.append('---------------------- SOLVER OPTIONS-----------------------------------------')
          self.finput.append('Option')
          self.finput.append('(-)')
          self.finput.append('help')
          self.finput.append(' integration_dt 0')
          self.finput.append(' kb_default 3.0e6')
          self.finput.append(' cb_default 3.0e5')
          self.finput.append(' wave_kinematics ')
          self.finput.append('inner_ftol 1e-5')
          self.finput.append('inner_gtol 1e-5')
          self.finput.append('inner_xtol 1e-5')
          self.finput.append('outer_tol 1e-3')
          self.finput.append(' pg_cooked 10000 1')
          self.finput.append(' outer_fd')
          self.finput.append(' outer_bd')
          self.finput.append(' outer_cd')
          self.finput.append(' inner_max_its 200')
          self.finput.append(' outer_max_its 600')
          # Repeat the details for the one mooring line multiple times
          angles = np.linspace(0, 360, n_connect+1)[1:-1]
          line = 'repeat'
          for degree in angles:
              line += (' %d' % degree)
          self.finput.append(line)
          self.finput.append(' krylov_accelerator 3')
          self.finput.append(' ref_position 0.0 0.0 0.0')
  
          
      def write_input_file(self, inputs, discrete_inputs):
          """Writes SOLVER OPTIONS section of input.map file, 
          which includes repeating node/line arrangement in even angular spacing around structure.
          
          INPUTS:
          ----------
          inputs        : dictionary of input parameters
          
          OUTPUTS  : none
          """
          # Unpack variables
          fairleadDepth = inputs['fairlead']
          R_fairlead    = inputs['fairlead_radius']
          R_anchor      = inputs['anchor_radius']
          n_connect     = int(inputs['number_of_mooring_connections'])
          n_lines       = int(inputs['mooring_lines_per_connection'])
          ntotal        = n_connect * n_lines
              
          
          # Open the map input file
          self.finput = []
  
          # Write the "Line Dictionary" section
          self.write_line_dictionary(inputs, discrete_inputs)
  
          # Write the "Node Properties" section
          self.write_node_properties_header()
          # One end on sea floor the other at fairlead
          self.write_node_properties(1, "VESSEL", R_fairlead, 0, -fairleadDepth)
          if n_lines > 1:
              angles  = np.linspace(0, 2*np.pi, ntotal+1)[:n_lines]
              angles -= np.mean(angles)
              anchorx = R_anchor * np.cos( angles )
              anchory = R_anchor * np.sin( angles )
              for k in range(n_lines):
                  self.write_node_properties(k+2, "FIX", anchorx[k], anchory[k], None)
          else:
              self.write_node_properties(2, "FIX", R_anchor, 0, None)
  
          # Write the "Line Properties" section
          for k in range(n_lines):
              self.write_line_properties(inputs, discrete_inputs, line_number=k+1, anchor_node=k+2, fairlead_node=1)
  
          # Write the "Solve Options" section
          self.write_solver_options(inputs)
          
          
      def runMAP(self, inputs, discrete_inputs, outputs):
          """Writes MAP input file, executes, and then queries MAP to find 
          maximum loading and displacement from vessel displacement around all 360 degrees
          
          INPUTS:
          ----------
          inputs   : dictionary of input parameters
          outputs : dictionary of output parameters
          
          OUTPUTS  : none (multiple unknown dictionary values set)
          """
          # Unpack variables
          rhoWater      = inputs['water_density']
          waterDepth    = inputs['water_depth']
          fairleadDepth = inputs['fairlead']
          Dmooring      = inputs['mooring_diameter']
          offset        = inputs['max_offset']
          heel          = inputs['operational_heel']
          gamma         = inputs['gamma_f']
          n_connect     = int(inputs['number_of_mooring_connections'])
          n_lines       = int(inputs['mooring_lines_per_connection'])
          ntotal        = n_connect * n_lines
  
          # Write the mooring system input file for this design
          self.write_input_file(inputs, discrete_inputs)
  
          # Initiate MAP++ for this design
          mymap = pyMAP( )
          #mymap.ierr = 0
          mymap.map_set_sea_depth(waterDepth)
          mymap.map_set_gravity(gravity)
          mymap.map_set_sea_density(rhoWater)
          mymap.read_list_input(self.finput)
          mymap.init( )
  
          # Get the stiffness matrix at neutral position
          mymap.displace_vessel(0, 0, 0, 0, 0, 0)
          mymap.update_states(0.0, 0)
          K = mymap.linear(1e-4) # Input finite difference epsilon
          outputs['mooring_stiffness'] = np.array( K )
          mymap.displace_vessel(0, 0, 0, 0, 0, 0)
          mymap.update_states(0.0, 0)
          
          # Get the vertical load on the structure and plotting data
          F_neutral = np.zeros((NLINES_MAX, 3))
          plotMat   = np.zeros((NLINES_MAX, NPTS_PLOT, 3))
          nptsMOI   = 100
          xyzpts    = np.zeros((ntotal, nptsMOI, 3)) # For MOI calculation
          for k in range(ntotal):
              (F_neutral[k,0], F_neutral[k,1], F_neutral[k,2]) = mymap.get_fairlead_force_3d(k)
              plotMat[k,:,0] = mymap.plot_x(k, NPTS_PLOT)
              plotMat[k,:,1] = mymap.plot_y(k, NPTS_PLOT)
              plotMat[k,:,2] = mymap.plot_z(k, NPTS_PLOT)
              xyzpts[k,:,0]  = mymap.plot_x(k, nptsMOI)
              xyzpts[k,:,1]  = mymap.plot_y(k, nptsMOI)
              xyzpts[k,:,2]  = mymap.plot_z(k, nptsMOI)
              if self.tlpFlag:
                  # Seems to be a bug in the plot arrays from MAP++ for plotting output with taut lines
                  plotMat[k,:,2] = np.linspace(-fairleadDepth, -waterDepth, NPTS_PLOT)
                  xyzpts[k,:,2]  = np.linspace(-fairleadDepth, -waterDepth, nptsMOI)
          outputs['mooring_neutral_load'] = F_neutral
          outputs['mooring_plot_matrix']  = plotMat
          
          # Fine line segment length, ds = sqrt(dx^2 + dy^2 + dz^2)
          xyzpts_dx = np.gradient(xyzpts[:,:,0], axis=1)
          xyzpts_dy = np.gradient(xyzpts[:,:,1], axis=1)
          xyzpts_dz = np.gradient(xyzpts[:,:,2], axis=1)
          xyzpts_ds = np.sqrt(xyzpts_dx**2 + xyzpts_dy**2 + xyzpts_dz**2)
  
          # Initialize inertia tensor integrands in https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor
          # Taking MOI relative to body centerline at fairlead depth
          r0 = np.array([0.0, 0.0, -fairleadDepth])
          R  = np.zeros((ntotal, nptsMOI, 6))
          for ii in range(nptsMOI):
              for k in range(ntotal):
                  r = xyzpts[k,ii,:] - r0
                  R[k,ii,:] = unassembleI(np.dot(r,r)*np.eye(3) - np.outer(r,r))
          Imat = self.wet_mass_per_length * np.trapz(R, x=xyzpts_ds[:,:,np.newaxis], axis=1)
          outputs['mooring_moments_of_inertia'] = np.abs( Imat.sum(axis=0) )
  
          # Get the restoring moment at maximum angle of heel
          # Since we don't know the substucture CG, have to just get the forces of the lines now and do the cross product later
          # We also want to allow for arbitraty wind direction and yaw of rotor relative to mooring lines, so we will compare
          # pitch and roll forces as extremes
          # TODO: This still isgn't quite the same as clocking the mooring lines in different directions,
          # which is what we want to do, but that requires multiple input files and solutions
          Fh = np.zeros((NLINES_MAX,3))
          mymap.displace_vessel(0, 0, 0, 0, heel, 0)
          mymap.update_states(0.0, 0)
          for k in range(ntotal):
              Fh[k][0], Fh[k][1], Fh[k][2] = mymap.get_fairlead_force_3d(k)
  
          outputs['operational_heel_restoring_force'] = Fh
  
          # Get angles by which to find the weakest line
          dangle  = 2.0
          angles  = np.deg2rad( np.arange(0.0, 360.0, dangle) )
          nangles = len(angles)
  
          # Get restoring force at weakest line at maximum allowable offset
          # Will global minimum always be along mooring angle?
          max_tension   = 0.0
          max_angle     = None
          min_angle     = None
          F_max_tension = None
          F_min         = np.inf
          T = np.zeros((NLINES_MAX,))
          F = np.zeros((NLINES_MAX,))
          # Loop around all angles to find weakest point
          for a in angles:
              # Unit vector and offset in x-y components
              idir  = np.array([np.cos(a), np.sin(a)])
              surge = offset * idir[0]
              sway  = offset * idir[1]
  
              # Get restoring force of offset at this angle
              mymap.displace_vessel(surge, sway, 0, 0, 0, 0) # 0s for z, angles
              mymap.update_states(0.0, 0)
              for k in range(ntotal):
                  # Force in x-y-z coordinates
                  fx,fy,fz = mymap.get_fairlead_force_3d(k)
                  T[k]     = np.sqrt(fx*fx + fy*fy + fz*fz)
                  # Total restoring force
                  F[k]     = np.dot([fx, fy], idir)
                  
              # Check if this is the weakest direction (highest tension)
              tempMax = T.max()
              if tempMax > max_tension:
                  max_tension   = tempMax
                  F_max_tension = F.sum()
                  max_angle     = a
              if F.sum() < F_min:
                  F_min     = F.sum()
                  min_angle = a
                  
          # Store the weakest restoring force when the vessel is offset the maximum amount
          outputs['max_offset_restoring_force'] = F_min
  
          # Check for good convergence
          if (plotMat[0,-1,-1] + fairleadDepth) > 1.0:
              outputs['axial_unity'] = 1e30
          else:
              outputs['axial_unity'] = gamma * max_tension / self.min_break_load
  
          mymap.end()
  
          
      def compute_cost(self, inputs, discrete_inputs, outputs):
          """Computes cost, based on mass scaling, of mooring system.
          
          INPUTS:
          ----------
          inputs   : dictionary of input parameters
          outputs : dictionary of output parameters
          
          OUTPUTS  : none (mooring_cost/mass unknown dictionary values set)
          """
          # Unpack variables
          rhoWater      = inputs['water_density']
          L_mooring     = inputs['mooring_line_length']
          anchorType    = discrete_inputs['anchor_type']
          costFact      = inputs['mooring_cost_factor']
          n_connect     = int(inputs['number_of_mooring_connections'])
          n_lines       = int(inputs['mooring_lines_per_connection'])
          ntotal        = n_connect * n_lines
          
          # Cost of anchors
          if type(anchorType) == type(''): anchorType = Anchor[anchorType.upper()]
          if anchorType == Anchor['DRAGEMBEDMENT']:
              anchor_rate = 1e-3 * self.min_break_load / gravity / 20*2000
          elif anchorType  == Anchor['SUCTIONPILE']:
              anchor_rate = 150000.* np.sqrt(1e-3*self.min_break_load/gravity/1250.)
          else:
              raise ValueError('Anchor Type must be DRAGEMBEDMENT or SUCTIONPILE')
          anchor_total = anchor_rate*ntotal
  
          # Cost of all of the mooring lines
          legs_total = ntotal * self.cost_per_length * L_mooring
  
          # Total summations
          outputs['anchor_cost']  = anchor_total
          outputs['mooring_cost'] = costFact*(legs_total + anchor_total)
          outputs['mooring_mass'] = self.wet_mass_per_length*L_mooring*ntotal
          outputs['number_of_mooring_lines'] = ntotal
  