diff --git a/docs/src/literate/C_customizablePreprocessing.jl b/docs/src/literate/C_customizablePreprocessing.jl
index 36cae866..05aa24fa 100644
--- a/docs/src/literate/C_customizablePreprocessing.jl
+++ b/docs/src/literate/C_customizablePreprocessing.jl
@@ -22,7 +22,7 @@ import PyPlot
 import OWENSOpenFASTWrappers
  
 
-path = runpath = "/home/runner/work/OWENS.jl/OWENS.jl/docs/src/literate" #splitdir(@__FILE__)[1]
+path = runpath = "/home/runner/work/OWENS.jl/OWENS.jl/docs/src/literate" #bld_height_numadsplitdir(@__FILE__)[1]
 
 Inp = OWENS.MasterInput("$runpath/sampleOWENS.yml")
 
@@ -408,11 +408,11 @@ if !AD15On
     delta_xs = shapeX[2:end] - shapeX[1:end-1]
     delta_zs = shapeZ[2:end] - shapeZ[1:end-1]
     delta3D = atan.(delta_xs./delta_zs)
-    delta3D_spl = FLOWMath.akima(shapeZ[1:end-1]./maximum(shapeZ[1:end-1]), delta3D,LinRange(0,1,length(numadIn_bld.span)-1))
+    delta3D_spl = OWENS.safeakima(shapeZ[1:end-1]./maximum(shapeZ[1:end-1]), delta3D,LinRange(0,1,length(numadIn_bld.span)-1))
     
     bld_height_numad = cumsum(diff(numadIn_bld.span).*(1.0.-abs.(sin.(delta3D_spl)))) # now convert the numad span to a height
 
-    chord = FLOWMath.akima(bld_height_numad./maximum(bld_height_numad), numadIn_bld.chord,LinRange(0,1,Nslices)) # now we can use it to access the numad data 
+    chord = OWENS.safeakima(bld_height_numad./maximum(bld_height_numad), numadIn_bld.chord,LinRange(minimum(bld_height_numad),1,Nslices)) # now we can use it to access the numad data 
     airfoils = fill("nothing",Nslices)
 
     for (iheight_numad,height_numad) in enumerate(bld_height_numad./maximum(bld_height_numad)) # Discretely assign the airfoils
@@ -453,7 +453,7 @@ if AD15On
 
     NumADBldNds = NumADStrutNds = 10 
 
-    bldchord_spl = FLOWMath.akima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,NumADBldNds))
+    bldchord_spl = OWENS.safeakima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,NumADBldNds))
 
     airfoil_filenames = fill("nothing",NumADBldNds) # Discretely assign the airfoils #TODO: separate out struts
     for (ispan_numad,span_numad) in enumerate(numadIn_bld.span./maximum(numadIn_bld.span))
@@ -470,7 +470,7 @@ if AD15On
         blade_Nnodes = [NumADBldNds for i=1:Nbld]
     else
         blade_filenames = [[blade_filename for i=1:Nbld];[lower_strut_filename for i=1:Nbld];[upper_strut_filename for i=1:Nbld]]
-        strutchord_spl = FLOWMath.akima(numadIn_strut.span./maximum(numadIn_strut.span), numadIn_strut.chord,LinRange(0,1,NumADStrutNds))
+        strutchord_spl = OWENS.safeakima(numadIn_strut.span./maximum(numadIn_strut.span), numadIn_strut.chord,LinRange(0,1,NumADStrutNds))
         blade_chords = [[bldchord_spl for i=1:Nbld];[strutchord_spl for i=1:Nbld];[strutchord_spl for i=1:Nbld]]
         blade_Nnodes = [[NumADBldNds for i=1:Nbld];[NumADStrutNds for i=1:Nbld];[NumADStrutNds for i=1:Nbld]]
     end
@@ -503,7 +503,7 @@ if AD15On
         ymesh = mymesh.y[strt_idx:end_idx]
         ADshapeX = sqrt.(xmesh.^2 .+ ymesh.^2)
         ADshapeX .-= ADshapeX[1] #get it starting at zero #TODO: make robust for blades that don't start at 0
-        ADshapeXspl = FLOWMath.akima(LinRange(0,H,length(ADshapeX)),ADshapeX,ADshapeZ)
+        ADshapeXspl = OWENS.safeakima(LinRange(0,H,length(ADshapeX)),ADshapeX,ADshapeZ)
         
         if iADBody<=Nbld #Note that the blades can be curved and are assumed to be oriented vertically
             BlSpn=ADshapeZ
diff --git a/examples/Optimization/optimizationExample.jl b/examples/Optimization/optimizationExample.jl
index cef0dc10..e60f01dd 100644
--- a/examples/Optimization/optimizationExample.jl
+++ b/examples/Optimization/optimizationExample.jl
@@ -90,7 +90,7 @@ function messyoptfun!(constraints,Vars)
     x_control_pts_grid = [0.0,0.25,0.5,0.75,1.0]
 
     x_grid = windio[:components][:blade][:internal_structure_2d_fem][:reference_axis][:x][:grid]
-    x_values = FLOWMath.akima(x_control_pts_grid,x_control_pts,x_grid)
+    x_values = safeakima(x_control_pts_grid,x_control_pts,x_grid)
 
     windio[:components][:blade][:internal_structure_2d_fem][:reference_axis][:x][:values] = x_values 
 
diff --git a/examples/SNL34m/34mSetup.jl b/examples/SNL34m/34mSetup.jl
index f6dea130..d356649c 100644
--- a/examples/SNL34m/34mSetup.jl
+++ b/examples/SNL34m/34mSetup.jl
@@ -17,7 +17,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 # Analytical
 # z_shape = [0.0, 0.4027099927689326, 0.8054199855378652, 1.2081299783067978, 1.6108399710757304, 2.013549963844663, 2.4162599566135956, 2.818969949382528, 3.221679942151461, 3.624389934920394, 4.027099927689326, 4.429809920458259, 4.832519913227191, 5.235229905996124, 5.637939898765056, 6.221683770581164, 6.821719178571302, 7.437583162221221, 8.068800548394838, 8.714884317956601, 9.375335981533686, 10.049645964128134, 10.737293998282153, 11.437749525493246, 12.305529745531164, 13.201631116890074, 14.122956691386433, 15.066322345375763, 16.028467784161606, 17.00606780965343, 17.995743812332336, 18.994075447807884, 19.997612457611687, 21.002886593374797, 22.006423603178604, 23.004755238654155, 23.99443124133306, 24.97203126682489, 25.934176705610724, 26.877542359600053, 27.79886793409642, 28.694969305455324, 29.562749525493246, 30.263205052704336, 30.950853086858352, 31.62516306945281, 32.285614733029895, 32.93169850259165, 33.56291588876527, 34.1787798724152, 34.77881528040533, 35.362559152221436, 35.821422444858186, 36.280285737494935, 36.739149030131685, 37.198012322768435, 37.656875615405184, 38.11573890804193, 38.57460220067868, 39.03346549331543, 39.49232878595218, 39.951192078588925, 40.410055371225674, 40.868918663862424, 41.32778195649917, 41.78664524913592]
@@ -102,8 +102,8 @@ end
 # PyPlot.figure()
 # PyPlot.plot(bld_OWENSPreCompinput[1].ynode,bld_OWENSPreCompinput[1].ynode)
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_OWENSPreComp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_OWENSPreComp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_OWENSPreComp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_OWENSPreComp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 
@@ -363,12 +363,12 @@ function runowens(model,feamodel,mymesh,myel,aeroForcesDMS,deformTurb;steady=tru
         # samplepts = numadIn_bld.span./maximum(numadIn_bld.span) #normalize #TODO: this is spanwise, while everything else is vertical-wise
         for its = 1:N_ts
             #TODO: there are strain values at each quad point, should be better than just choosing one
-            eps_x[ibld,its,:] = epsilon_x_hist[1,start:stop,its]#FLOWMath.akima(x,epsilon_x_hist[1,start:stop,its],samplepts)
-            eps_z[ibld,its,:] = meanepsilon_z_hist[1,start:stop,its]#FLOWMath.akima(x,meanepsilon_z_hist[1,start:stop,its],samplepts)
-            eps_y[ibld,its,:] = meanepsilon_y_hist[1,start:stop,its]#FLOWMath.akima(x,meanepsilon_y_hist[1,start:stop,its],samplepts)
-            kappa_x[ibld,its,:] = kappa_x_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_x_hist[1,start:stop,its],samplepts)
-            kappa_y[ibld,its,:] = kappa_y_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_y_hist[1,start:stop,its],samplepts)
-            kappa_z[ibld,its,:] = kappa_z_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_z_hist[1,start:stop,its],samplepts)
+            eps_x[ibld,its,:] = epsilon_x_hist[1,start:stop,its]#safeakima(x,epsilon_x_hist[1,start:stop,its],samplepts)
+            eps_z[ibld,its,:] = meanepsilon_z_hist[1,start:stop,its]#safeakima(x,meanepsilon_z_hist[1,start:stop,its],samplepts)
+            eps_y[ibld,its,:] = meanepsilon_y_hist[1,start:stop,its]#safeakima(x,meanepsilon_y_hist[1,start:stop,its],samplepts)
+            kappa_x[ibld,its,:] = kappa_x_hist[1,start:stop,its]#safeakima(x,kappa_x_hist[1,start:stop,its],samplepts)
+            kappa_y[ibld,its,:] = kappa_y_hist[1,start:stop,its]#safeakima(x,kappa_y_hist[1,start:stop,its],samplepts)
+            kappa_z[ibld,its,:] = kappa_z_hist[1,start:stop,its]#safeakima(x,kappa_z_hist[1,start:stop,its],samplepts)
         end
     end
 
diff --git a/examples/SNL34m/Example34mHVAWTConfiguration.jl b/examples/SNL34m/Example34mHVAWTConfiguration.jl
index 9e2eb244..7b36db74 100644
--- a/examples/SNL34m/Example34mHVAWTConfiguration.jl
+++ b/examples/SNL34m/Example34mHVAWTConfiguration.jl
@@ -81,9 +81,9 @@ SNL34m_5_3_Torque = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Data/5.3_T
 
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 offsetTime = 20.0 # seconds
 tocp = [0.0;new_t.+offsetTime; 1e6]
@@ -99,7 +99,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -234,8 +234,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 
 flatwise_stress1grav = (kappa_y_grav[1,end-1,2:end].* thickness .+ 0*eps_x_grav[1,end-1,2:end]) .* Ealuminum
 flatwise_stress2grav = (kappa_y_grav[2,end-1,1:end-1].* thickness .+ 0*eps_x_grav[2,end-1,1:end-1]) .* Ealuminum
diff --git a/examples/SNL34m/Fig4.1Gravity_Fig4.4_Steady40RPMCentrifugalGX.jl b/examples/SNL34m/Fig4.1Gravity_Fig4.4_Steady40RPMCentrifugalGX.jl
index b06e1a81..c7bc0e86 100644
--- a/examples/SNL34m/Fig4.1Gravity_Fig4.4_Steady40RPMCentrifugalGX.jl
+++ b/examples/SNL34m/Fig4.1Gravity_Fig4.4_Steady40RPMCentrifugalGX.jl
@@ -57,7 +57,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -130,8 +130,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 PyPlot.figure()
diff --git a/examples/SNL34m/Fig5_11_EmergStop_torque_flatwise.jl b/examples/SNL34m/Fig5_11_EmergStop_torque_flatwise.jl
index 125ea498..7fb2291a 100644
--- a/examples/SNL34m/Fig5_11_EmergStop_torque_flatwise.jl
+++ b/examples/SNL34m/Fig5_11_EmergStop_torque_flatwise.jl
@@ -98,12 +98,12 @@ SNL34m_5_11_Vinf = SNL34m_5_11_Vinf[1:2:end,:]
 SNL34m_5_11_RPM[SNL34m_5_11_RPM[:,2].<0.0,2] .*= 0
 
 new_t = LinRange(SNL34m_5_11_RPM[1,1],SNL34m_5_11_RPM[end,1],100).-SNL34m_5_11_RPM[1,1].+1e-6
-new_RPM = FLOWMath.akima(SNL34m_5_11_RPM[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_11_RPM[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_RPM[:,2],new_t)
 
-new_Torque = FLOWMath.akima(SNL34m_5_11_Torque[:,1],SNL34m_5_11_Torque[:,2],new_t)
+new_Torque = safeakima(SNL34m_5_11_Torque[:,1],SNL34m_5_11_Torque[:,2],new_t)
 
 t_Vinf = LinRange(SNL34m_5_11_Vinf[1,1],SNL34m_5_11_Vinf[end,1],60).-SNL34m_5_11_RPM[1,1].+1e-6
-Vinf_spec = FLOWMath.akima(SNL34m_5_11_Vinf[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_Vinf[:,2],t_Vinf)
+Vinf_spec = safeakima(SNL34m_5_11_Vinf[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_Vinf[:,2],t_Vinf)
 
 offsetTime = 20.0
 t_Vinf = [0;t_Vinf.+offsetTime;t_Vinf[end]*10]
@@ -127,7 +127,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -200,8 +200,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 PyPlot.figure()
diff --git a/examples/SNL34m/Fig5_4_Fig5_3_normaloperation_torque_flatwise_generatorControl.jl b/examples/SNL34m/Fig5_4_Fig5_3_normaloperation_torque_flatwise_generatorControl.jl
index 285cde56..36dee068 100644
--- a/examples/SNL34m/Fig5_4_Fig5_3_normaloperation_torque_flatwise_generatorControl.jl
+++ b/examples/SNL34m/Fig5_4_Fig5_3_normaloperation_torque_flatwise_generatorControl.jl
@@ -78,11 +78,11 @@ PyPlot.xlim([0,100.0])
 # PyPlot.savefig("$(path)/../figs/NormalOperation_Vinf.pdf",transparent = true)
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-new_Torque = FLOWMath.akima(SNL34m_5_3_Torque[:,1],SNL34m_5_3_Torque[:,2],new_t)
+new_Torque = safeakima(SNL34m_5_3_Torque[:,1],SNL34m_5_3_Torque[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 t_Vinf = [0;new_t;1e6]
 Vinf_spec = [Vinf_spec[1];Vinf_spec;Vinf_spec[end]]
@@ -121,7 +121,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -194,8 +194,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 
diff --git a/examples/SNL34m/SNL34mVAWTNormalOperation.jl b/examples/SNL34m/SNL34mVAWTNormalOperation.jl
index ba7fd385..1dca2737 100644
--- a/examples/SNL34m/SNL34mVAWTNormalOperation.jl
+++ b/examples/SNL34m/SNL34mVAWTNormalOperation.jl
@@ -81,9 +81,9 @@ SNL34m_5_3_Torque = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Data/5.3_T
 
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = OWENS.safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = OWENS.safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 offsetTime = 20.0 # seconds
 tocp = [0.0;new_t.+offsetTime; 1e6]
@@ -99,7 +99,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = OWENS.safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -130,6 +130,7 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections,AD15bldNdIdxRng,
     R = Blade_Radius,
     shapeZ,
     shapeX,
+    shapeY = zero(shapeX),
     ifw,
     delta_t,
     numTS,
@@ -353,8 +354,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = OWENS.safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = OWENS.safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 
 flatwise_stress1grav = (kappa_y_grav[1,end-1,2:end].* thickness .+ 0*eps_x_grav[1,end-1,2:end]) .* Ealuminum
 flatwise_stress2grav = (kappa_y_grav[2,end-1,1:end-1].* thickness .+ 0*eps_x_grav[2,end-1,1:end-1]) .* Ealuminum
diff --git a/examples/SNL34m/SNL34m_Inputs.yml b/examples/SNL34m/SNL34m_Inputs.yml
index dd1d59c1..925b210a 100644
--- a/examples/SNL34m/SNL34m_Inputs.yml
+++ b/examples/SNL34m/SNL34m_Inputs.yml
@@ -23,7 +23,7 @@ turbulentInflow:
 controlParameters:
   controlStrategy: constantRPM # TODO: incorporate the others
   RPM: 34.0 #RPM
-  numTS: 10 #
+  numTS: 5000 #
   delta_t: 0.05 # s
 
 
@@ -35,9 +35,9 @@ AeroParameters:
   adi_rootname: "/SNL34m"
   
 structuralParameters:
-  structuralModel: ROM #GX, TNB, ROM
+  structuralModel: GX #GX, TNB, ROM
   nonlinear: false #TODO: propogate
-  ntelem: 10 #tower elements in each 
+  ntelem: 20 #tower elements in each 
   nbelem: 60 #blade elements in each 
   ncelem: 10 #central cable elements in each if turbineType is ARCUS
   nselem: 5 #strut elements in each if turbineType has struts
diff --git a/examples/SNL34m/TSR_CP_2way_Aerostructural.jl b/examples/SNL34m/TSR_CP_2way_Aerostructural.jl
index 2879db04..559efa54 100644
--- a/examples/SNL34m/TSR_CP_2way_Aerostructural.jl
+++ b/examples/SNL34m/TSR_CP_2way_Aerostructural.jl
@@ -75,9 +75,9 @@ SNL34m_5_3_Torque = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Data/5.3_T
 
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 offsetTime = 20.0 # seconds
 tocp = [0.0;new_t.+offsetTime; 1e6]
@@ -93,7 +93,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
diff --git a/examples/dev/OptimizationExample.jl b/examples/dev/OptimizationExample.jl
index 9a40f6af..4978b05f 100644
--- a/examples/dev/OptimizationExample.jl
+++ b/examples/dev/OptimizationExample.jl
@@ -107,7 +107,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
diff --git a/examples/dev/RM2.jl b/examples/dev/RM2.jl
index 086525fa..4ec7f9cd 100644
--- a/examples/dev/RM2.jl
+++ b/examples/dev/RM2.jl
@@ -104,7 +104,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
diff --git a/src/AeroMapping.jl b/src/AeroMapping.jl
index b6906fa2..c6d8b43b 100644
--- a/src/AeroMapping.jl
+++ b/src/AeroMapping.jl
@@ -115,12 +115,12 @@ function mapACDMS(t,azi_j,mesh,el,advanceTurb;numAeroTS = 1,alwaysrecalc=true,ou
     for i=1:NBlade
         for j=1:numAeroTS
             # AC and DMS calculate inbetween aero slices, so we add the 0 and 1 normed values here to ensure we don't extrapolate
-            struct_N[i,j,:] = FLOWMath.akima([0.0;spanLocNorm[i,:];1.0],[0.0;N[i,j,:];0.0],structuralSpanLocNorm[i,:])
-            struct_T[i,j,:] = FLOWMath.akima([0.0;spanLocNorm[i,:];1.0],[0.0;T[i,j,:];0.0],structuralSpanLocNorm[i,:])
-            struct_M25[i,j,:] = FLOWMath.akima([0.0;spanLocNorm[i,:];1.0],[0.0;M25[i,j,:];0.0],structuralSpanLocNorm[i,:])
-            struct_X[i,j,:] = FLOWMath.akima([0.0;spanLocNorm[i,:];1.0],[0.0;X[i,j,:];0.0],structuralSpanLocNorm[i,:])
-            struct_Y[i,j,:] = FLOWMath.akima([0.0;spanLocNorm[i,:];1.0],[0.0;Y[i,j,:];0.0],structuralSpanLocNorm[i,:])
-            struct_Z[i,j,:] = FLOWMath.akima([0.0;spanLocNorm[i,:];1.0],[0.0;Z[i,j,:];0.0],structuralSpanLocNorm[i,:])
+            struct_N[i,j,:] = safeakima([0.0;spanLocNorm[i,:];1.0],[0.0;N[i,j,:];0.0],structuralSpanLocNorm[i,:])
+            struct_T[i,j,:] = safeakima([0.0;spanLocNorm[i,:];1.0],[0.0;T[i,j,:];0.0],structuralSpanLocNorm[i,:])
+            struct_M25[i,j,:] = safeakima([0.0;spanLocNorm[i,:];1.0],[0.0;M25[i,j,:];0.0],structuralSpanLocNorm[i,:])
+            struct_X[i,j,:] = safeakima([0.0;spanLocNorm[i,:];1.0],[0.0;X[i,j,:];0.0],structuralSpanLocNorm[i,:])
+            struct_Y[i,j,:] = safeakima([0.0;spanLocNorm[i,:];1.0],[0.0;Y[i,j,:];0.0],structuralSpanLocNorm[i,:])
+            struct_Z[i,j,:] = safeakima([0.0;spanLocNorm[i,:];1.0],[0.0;Z[i,j,:];0.0],structuralSpanLocNorm[i,:])
         end
     end
 
diff --git a/src/PostProcessing.jl b/src/PostProcessing.jl
index 4e6daf35..2d043892 100644
--- a/src/PostProcessing.jl
+++ b/src/PostProcessing.jl
@@ -458,7 +458,7 @@ function calcSF(stress,SF_ult,SF_buck,lencomposites_span,plyprops,
                     reverse!(Log_SN_cycles2Fail1)
                 end
 
-                logcycles2fail = FLOWMath.akima(SN_stressMpa1.*1e6,Log_SN_cycles2Fail1,stress_levels)
+                logcycles2fail = safeakima(SN_stressMpa1.*1e6,Log_SN_cycles2Fail1,stress_levels)
                 cycles2fail = 10.0 .^ logcycles2fail
                 damage_layers[ilayer] = sum(cyclesatStressRanges./cycles2fail)
             end
@@ -695,17 +695,17 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
         end
 
         mesh_span_bld = mymesh.z[startN:stopN].-mymesh.z[startN]
-        mesh_span_bld_el = FLOWMath.akima(LinRange(0,1,length(mesh_span_bld)),mesh_span_bld,LinRange(0,1,stopE-startE+1))
+        mesh_span_bld_el = safeakima(LinRange(0,1,length(mesh_span_bld)),mesh_span_bld,LinRange(0,1,stopE-startE+1))
         composites_span_bld = numadIn_bld.span/maximum(numadIn_bld.span)*maximum(mesh_span_bld_el) #TODO: this assumes struts and blades are straight, should be mapped for all potential cases, also need to input hub offset
 
         for its = 1:N_ts
             # Interpolate to the composite inputs #TODO: verify node vs el in strain
-            eps_x_bld[ibld,its,:] = FLOWMath.akima(mesh_span_bld_el,epsilon_x_hist[1,startE:stopE,its],composites_span_bld)
-            eps_z_bld[ibld,its,:] = FLOWMath.akima(mesh_span_bld_el,meanepsilon_z_hist[1,startE:stopE,its],composites_span_bld)
-            eps_y_bld[ibld,its,:] = FLOWMath.akima(mesh_span_bld_el,meanepsilon_y_hist[1,startE:stopE,its],composites_span_bld)
-            kappa_x_bld[ibld,its,:] = FLOWMath.akima(mesh_span_bld_el,kappa_x_hist[1,startE:stopE,its],composites_span_bld)
-            kappa_y_bld[ibld,its,:] = FLOWMath.akima(mesh_span_bld_el,kappa_y_hist[1,startE:stopE,its],composites_span_bld)
-            kappa_z_bld[ibld,its,:] = FLOWMath.akima(mesh_span_bld_el,kappa_z_hist[1,startE:stopE,its],composites_span_bld)
+            eps_x_bld[ibld,its,:] = safeakima(mesh_span_bld_el,epsilon_x_hist[1,startE:stopE,its],composites_span_bld)
+            eps_z_bld[ibld,its,:] = safeakima(mesh_span_bld_el,meanepsilon_z_hist[1,startE:stopE,its],composites_span_bld)
+            eps_y_bld[ibld,its,:] = safeakima(mesh_span_bld_el,meanepsilon_y_hist[1,startE:stopE,its],composites_span_bld)
+            kappa_x_bld[ibld,its,:] = safeakima(mesh_span_bld_el,kappa_x_hist[1,startE:stopE,its],composites_span_bld)
+            kappa_y_bld[ibld,its,:] = safeakima(mesh_span_bld_el,kappa_y_hist[1,startE:stopE,its],composites_span_bld)
+            kappa_z_bld[ibld,its,:] = safeakima(mesh_span_bld_el,kappa_z_hist[1,startE:stopE,its],composites_span_bld)
         end
     end
 
@@ -742,17 +742,17 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
             end
 
             mesh_span_strut = abs.(mymesh.z[startN:stopN].-mymesh.z[startN])
-            mesh_span_strut_el = FLOWMath.akima(LinRange(0,1,length(mesh_span_strut)),mesh_span_strut,LinRange(0,1,stopE-startE+1))
+            mesh_span_strut_el = safeakima(LinRange(0,1,length(mesh_span_strut)),mesh_span_strut,LinRange(0,1,stopE-startE+1))
             composites_span_strut = numadIn_strut.span/maximum(numadIn_strut.span)*maximum(mesh_span_strut_el) #TODO: this assumes struts and blades are straight, should be mapped for all potential cases, also need to input hub offset
 
             for its = 1:N_ts
                 # Interpolate to the composite inputs #TODO: verify node vs el in strain
-                eps_x_strut[istrut,its,:] = FLOWMath.akima(mesh_span_strut_el,epsilon_x_hist[1,startE:stopE,its],composites_span_strut)
-                eps_z_strut[istrut,its,:] = FLOWMath.akima(mesh_span_strut_el,meanepsilon_z_hist[1,startE:stopE,its],composites_span_strut)
-                eps_y_strut[istrut,its,:] = FLOWMath.akima(mesh_span_strut_el,meanepsilon_y_hist[1,startE:stopE,its],composites_span_strut)
-                kappa_x_strut[istrut,its,:] = FLOWMath.akima(mesh_span_strut_el,kappa_x_hist[1,startE:stopE,its],composites_span_strut)
-                kappa_y_strut[istrut,its,:] = FLOWMath.akima(mesh_span_strut_el,kappa_y_hist[1,startE:stopE,its],composites_span_strut)
-                kappa_z_strut[istrut,its,:] = FLOWMath.akima(mesh_span_strut_el,kappa_z_hist[1,startE:stopE,its],composites_span_strut)
+                eps_x_strut[istrut,its,:] = safeakima(mesh_span_strut_el,epsilon_x_hist[1,startE:stopE,its],composites_span_strut)
+                eps_z_strut[istrut,its,:] = safeakima(mesh_span_strut_el,meanepsilon_z_hist[1,startE:stopE,its],composites_span_strut)
+                eps_y_strut[istrut,its,:] = safeakima(mesh_span_strut_el,meanepsilon_y_hist[1,startE:stopE,its],composites_span_strut)
+                kappa_x_strut[istrut,its,:] = safeakima(mesh_span_strut_el,kappa_x_hist[1,startE:stopE,its],composites_span_strut)
+                kappa_y_strut[istrut,its,:] = safeakima(mesh_span_strut_el,kappa_y_hist[1,startE:stopE,its],composites_span_strut)
+                kappa_z_strut[istrut,its,:] = safeakima(mesh_span_strut_el,kappa_z_hist[1,startE:stopE,its],composites_span_strut)
             end
         end
     end
@@ -775,15 +775,15 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
     x = x.-x[1] #zero
     x = x./x[end] #normalize
     mesh_span_twr = mymesh.z[start:stop].-mymesh.z[start]
-    composites_span_twr = FLOWMath.akima(LinRange(0,1,length(mesh_span_twr)),mesh_span_twr,LinRange(0,1,length(twr_precompinput)))
+    composites_span_twr = safeakima(LinRange(0,1,length(mesh_span_twr)),mesh_span_twr,LinRange(0,1,length(twr_precompinput)))
     for its = 1:N_ts
         # Interpolate to the composite inputs
-        eps_x_twr[1,its,:] = FLOWMath.akima(mesh_span_twr,epsilon_x_hist[1,start:stop,its],composites_span_twr)
-        eps_z_twr[1,its,:] = FLOWMath.akima(mesh_span_twr,meanepsilon_z_hist[1,start:stop,its],composites_span_twr)
-        eps_y_twr[1,its,:] = FLOWMath.akima(mesh_span_twr,meanepsilon_y_hist[1,start:stop,its],composites_span_twr)
-        kappa_x_twr[1,its,:] = FLOWMath.akima(mesh_span_twr,kappa_x_hist[1,start:stop,its],composites_span_twr)
-        kappa_y_twr[1,its,:] = FLOWMath.akima(mesh_span_twr,kappa_y_hist[1,start:stop,its],composites_span_twr)
-        kappa_z_twr[1,its,:] = FLOWMath.akima(mesh_span_twr,kappa_z_hist[1,start:stop,its],composites_span_twr)
+        eps_x_twr[1,its,:] = safeakima(mesh_span_twr,epsilon_x_hist[1,start:stop,its],composites_span_twr)
+        eps_z_twr[1,its,:] = safeakima(mesh_span_twr,meanepsilon_z_hist[1,start:stop,its],composites_span_twr)
+        eps_y_twr[1,its,:] = safeakima(mesh_span_twr,meanepsilon_y_hist[1,start:stop,its],composites_span_twr)
+        kappa_x_twr[1,its,:] = safeakima(mesh_span_twr,kappa_x_hist[1,start:stop,its],composites_span_twr)
+        kappa_y_twr[1,its,:] = safeakima(mesh_span_twr,kappa_y_hist[1,start:stop,its],composites_span_twr)
+        kappa_z_twr[1,its,:] = safeakima(mesh_span_twr,kappa_z_hist[1,start:stop,its],composites_span_twr)
 
     end
 
diff --git a/src/SetupTurbine.jl b/src/SetupTurbine.jl
index e77764f3..ee005f44 100644
--- a/src/SetupTurbine.jl
+++ b/src/SetupTurbine.jl
@@ -427,16 +427,16 @@ function setupOWENS(OWENSAero,path;
         delta_xs = shapeX[2:end] - shapeX[1:end-1]
         delta_zs = shapeZ[2:end] - shapeZ[1:end-1]
         delta3D = atan.(delta_xs./delta_zs)
-        delta3D_spl = FLOWMath.akima(shapeZ[1:end-1]./maximum(shapeZ[1:end-1]), delta3D,LinRange(0,1,length(numadIn_bld.span)-1))
+        delta3D_spl = safeakima(shapeZ[1:end-1]./maximum(shapeZ[1:end-1]), delta3D,LinRange(0,1,length(numadIn_bld.span)-1))
         # now convert the numad span to a height
         bld_height_numad = cumsum(diff(numadIn_bld.span).*(1.0.-abs.(sin.(delta3D_spl))))
-
+        bld_height_numad_unit = bld_height_numad./maximum(bld_height_numad)
         # now we can use it to access the numad data 
-        chord = FLOWMath.akima(bld_height_numad./maximum(bld_height_numad), numadIn_bld.chord,LinRange(0,1,Nslices))
+        chord = safeakima(bld_height_numad_unit, numadIn_bld.chord,LinRange(bld_height_numad_unit[1],1,Nslices))
         airfoils = fill("nothing",Nslices)
 
         # Discretely assign the airfoils
-        for (iheight_numad,height_numad) in enumerate(bld_height_numad./maximum(bld_height_numad))
+        for (iheight_numad,height_numad) in enumerate(bld_height_numad_unit)
             for (iheight,height_slices) in enumerate(collect(LinRange(0,1,Nslices)))
                 if airfoils[iheight]=="nothing" && height_slices<=height_numad
                     airfoils[iheight] = "$(numadIn_bld.airfoil[iheight_numad]).dat"
@@ -473,7 +473,7 @@ function setupOWENS(OWENSAero,path;
 
         NumADBldNds = NumADStrutNds = 10 
 
-        bldchord_spl = FLOWMath.akima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,NumADBldNds))
+        bldchord_spl = safeakima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,NumADBldNds))
 
         # Discretely assign the blade airfoils based on the next closest neighbor
         bld_airfoil_filenames = fill("nothing",NumADBldNds) #TODO: cable drag?
@@ -498,7 +498,7 @@ function setupOWENS(OWENSAero,path;
             airfoil_filenames = collect(Iterators.flatten([bld_airfoil_filenames for i=1:Nbld]))
             
             for istrut = 1:Nstrutperbld
-                strutchord_spl = FLOWMath.akima(numadIn_strut[istrut].span./maximum(numadIn_strut[istrut].span), numadIn_strut[istrut].chord,LinRange(0,1,NumADStrutNds))
+                strutchord_spl = safeakima(numadIn_strut[istrut].span./maximum(numadIn_strut[istrut].span), numadIn_strut[istrut].chord,LinRange(0,1,NumADStrutNds))
                 for ibld = 1:Nbld
                     blade_filenames = [blade_filenames;"$path/strut$(istrut)_bld$ibld.dat"]
                     blade_chords = [blade_chords;[strutchord_spl]]
@@ -549,7 +549,7 @@ function setupOWENS(OWENSAero,path;
             ymesh = mymesh.y[strt_idx:end_idx]
             ADshapeX = sqrt.(xmesh.^2 .+ ymesh.^2)
             ADshapeX .-= ADshapeX[1] #get it starting at zero #TODO: make robust for blades that don't start at 0
-            ADshapeXspl = FLOWMath.akima(LinRange(0,H,length(ADshapeX)),ADshapeX,ADshapeZ)
+            ADshapeXspl = safeakima(LinRange(0,H,length(ADshapeX)),ADshapeX,ADshapeZ)
             
             if iADBody<=Nbld #&& !bladefileissaved#Note that the blades can be curved and are assumed to be oriented vertically
                 # bladefileissaved = true
@@ -564,16 +564,16 @@ function setupOWENS(OWENSAero,path;
                 #TODO: reevalueate these equations and make sure they are robust against varying designs
                 BlCrvACinput = -ymesh.*sin(bladeangle).+xmesh.*cos(bladeangle)
                 BlCrvACinput = BlCrvACinput .- BlCrvACinput[1]
-                BlSwpAC = -FLOWMath.akima(LinRange(0,H,length(BlCrvACinput)),BlCrvACinput,ADshapeZ)
+                BlSwpAC = -safeakima(LinRange(0,H,length(BlCrvACinput)),BlCrvACinput,ADshapeZ)
 
                 BlSwpACinput = xmesh.*sin(bladeangle).+ymesh.*cos(bladeangle)
                 BlSwpACinput = BlSwpACinput .- BlSwpACinput[1]
-                BlCrvAC = FLOWMath.akima(LinRange(0,H,length(BlSwpACinput)),BlSwpACinput,ADshapeZ)
+                BlCrvAC = safeakima(LinRange(0,H,length(BlSwpACinput)),BlSwpACinput,ADshapeZ)
 
                 BlCrvAng = zeros(blade_Nnodes[iADBody])
 
                 BlTwistinput =(blade_twist.-blade_twist[1])*180/pi
-                BlTwist = FLOWMath.akima(LinRange(0,H,length(BlTwistinput)),BlTwistinput,ADshapeZ)
+                BlTwist = safeakima(LinRange(0,H,length(BlTwistinput)),BlTwistinput,ADshapeZ)
 
                 BlChord=blade_chords[iADBody]
 
@@ -736,9 +736,9 @@ function setupOWENShawt(OWENSAero,path;
     omega = RPM / 60 * 2 * pi
     tsr = omega*R/Vinf
 
-    shapeX_spline = FLOWMath.Akima(shapeZ, shapeX)
-    bladelen = sum(sqrt.((shapeX[2:end].-shapeX[1:end-1]).^2 .+ (shapeZ[2:end].-shapeZ[1:end-1]).^2 ))
-    # println("bladelen 161.06: $bladelen")
+    # shapeX_spline = FLOWMath.Akima(shapeZ, shapeX)
+    # bladelen = sum(sqrt.((shapeX[2:end].-shapeX[1:end-1]).^2 .+ (shapeZ[2:end].-shapeZ[1:end-1]).^2 ))
+    # # println("bladelen 161.06: $bladelen")
     h_frac = (shapeZ[2:end] - shapeZ[1:end-1])./shapeZ[end];
     h_elem = (shapeZ[2:end] - shapeZ[1:end-1])
     h = (shapeZ[2:end] + shapeZ[1:end-1])/2.0;
@@ -933,7 +933,7 @@ function setupOWENShawt(OWENSAero,path;
     ### Set up aero forces
     #########################################
     # println("Initialize Aerodynamics")
-    # chord_spl = FLOWMath.akima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,Nslices))
+    # chord_spl = safeakima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,Nslices))
     # OWENSAero.setupTurb(shapeX,shapeZ,B,chord_spl,tsr,Vinf;AModel,DSModel,
     # afname = "$path/Airfoils/NACA_0021.dat", #TODO: map to the numad input
     # ifw,
diff --git a/src/Unsteady.jl b/src/Unsteady.jl
index fc212499..48feb862 100644
--- a/src/Unsteady.jl
+++ b/src/Unsteady.jl
@@ -290,7 +290,7 @@ function Unsteady(inputs;topModel=nothing,topMesh=nothing,topEl=nothing,
                 if inputs.generatorOn
                         if inputs.useGeneratorFunction
                             specifiedOmega,_,_ = omegaSpecCheck(t[i]+delta_t,inputs.tocp,inputs.Omegaocp,delta_t)
-                            newVinf = FLOWMath.akima(inputs.tocp_Vinf,inputs.Vinfocp,t[i])
+                            newVinf = safeakima(inputs.tocp_Vinf,inputs.Vinfocp,t[i])
                             if isnothing(userDefinedGenerator)
                                 genTorqueHSS0,integrator_j,controlnamecurrent = internaluserDefinedGenerator(newVinf,t[i],azi_j,Omega_j,OmegaHist[i],OmegaDot_j,OmegaDotHist[i],delta_t,integrator,specifiedOmega) #;operPhase
                             else
@@ -938,7 +938,7 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
     if inputs.tocp_Vinf == -1
         newVinf = -1
     else
-        newVinf = FLOWMath.akima(inputs.tocp_Vinf,inputs.Vinfocp,t_i)
+        newVinf = safeakima(inputs.tocp_Vinf,inputs.Vinfocp,t_i)
     end
 
     if inputs.aeroLoadsOn > 1
diff --git a/src/Unsteady_Land.jl b/src/Unsteady_Land.jl
index 39d4349b..9d286e59 100644
--- a/src/Unsteady_Land.jl
+++ b/src/Unsteady_Land.jl
@@ -303,7 +303,7 @@ function Unsteady_Land(inputs;topModel=nothing,topMesh=nothing,topEl=nothing,
             if inputs.generatorOn
                     if inputs.useGeneratorFunction
                         specifiedOmega,_,_ = omegaSpecCheck(t[i]+topdata.delta_t,inputs.tocp,inputs.Omegaocp,topdata.delta_t)
-                        newVinf = FLOWMath.akima(inputs.tocp_Vinf,inputs.Vinfocp,t[i]) #TODO: ifw sampling of same file as aerodyn
+                        newVinf = safeakima(inputs.tocp_Vinf,inputs.Vinfocp,t[i]) #TODO: ifw sampling of same file as aerodyn
                         if isnothing(userDefinedGenerator)
                             genTorqueHSS0,topdata.integrator_j,controlnamecurrent = internaluserDefinedGenerator(newVinf,t[i],topdata.azi_j,topdata.Omega_j,topdata.OmegaHist[i],topdata.OmegaDot_j,topdata.OmegaDotHist[i],topdata.delta_t,topdata.integrator,specifiedOmega) #;operPhase
                         else
@@ -742,12 +742,12 @@ function run34m(inputs,feamodel,mymesh,myel,aeroForces,deformAero;steady=true,sy
         # samplepts = numadIn_bld.span./maximum(numadIn_bld.span) #normalize #TODO: this is spanwise, while everything else is vertical-wise
         for its = 1:N_ts
             #TODO: there are strain values at each quad point, should be better than just choosing one
-            eps_x[ibld,its,:] = epsilon_x_hist[1,start:stop,its]#FLOWMath.akima(x,epsilon_x_hist[1,start:stop,its],samplepts)
-            eps_z[ibld,its,:] = meanepsilon_z_hist[1,start:stop,its]#FLOWMath.akima(x,meanepsilon_z_hist[1,start:stop,its],samplepts)
-            eps_y[ibld,its,:] = meanepsilon_y_hist[1,start:stop,its]#FLOWMath.akima(x,meanepsilon_y_hist[1,start:stop,its],samplepts)
-            kappa_x[ibld,its,:] = kappa_x_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_x_hist[1,start:stop,its],samplepts)
-            kappa_y[ibld,its,:] = kappa_y_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_y_hist[1,start:stop,its],samplepts)
-            kappa_z[ibld,its,:] = kappa_z_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_z_hist[1,start:stop,its],samplepts)
+            eps_x[ibld,its,:] = epsilon_x_hist[1,start:stop,its]#safeakima(x,epsilon_x_hist[1,start:stop,its],samplepts)
+            eps_z[ibld,its,:] = meanepsilon_z_hist[1,start:stop,its]#safeakima(x,meanepsilon_z_hist[1,start:stop,its],samplepts)
+            eps_y[ibld,its,:] = meanepsilon_y_hist[1,start:stop,its]#safeakima(x,meanepsilon_y_hist[1,start:stop,its],samplepts)
+            kappa_x[ibld,its,:] = kappa_x_hist[1,start:stop,its]#safeakima(x,kappa_x_hist[1,start:stop,its],samplepts)
+            kappa_y[ibld,its,:] = kappa_y_hist[1,start:stop,its]#safeakima(x,kappa_y_hist[1,start:stop,its],samplepts)
+            kappa_z[ibld,its,:] = kappa_z_hist[1,start:stop,its]#safeakima(x,kappa_z_hist[1,start:stop,its],samplepts)
         end
     end
 
diff --git a/src/fileio.jl b/src/fileio.jl
index 216f2aeb..3d08bfb9 100644
--- a/src/fileio.jl
+++ b/src/fileio.jl
@@ -42,7 +42,7 @@ function readNuMadGeomCSV(NuMad_geom_file::OrderedCollections.OrderedDict{Symbol
     airfoil = Array{String,1}(undef,length(span))
     airfoil_names = sec_Dict[:outer_shape_bem][:airfoil_position][:labels]
     # spline the airfoils used to the current overall grid, then round to enable mapping to the discrete airfoil inputs
-    airfoil_station_numbers = round.(Int,FLOWMath.akima(airfoil_grid,collect(1:length(airfoil_names)),span))
+    airfoil_station_numbers = round.(Int,safeakima(airfoil_grid,collect(1:length(airfoil_names)),span))
     # Now map the resulting airfoils to the new grid
     for istation = 1:length(airfoil_station_numbers)
         for iaf = 1:length(airfoil_names)
@@ -65,15 +65,15 @@ function readNuMadGeomCSV(NuMad_geom_file::OrderedCollections.OrderedDict{Symbol
 
     twist_grid = sec_Dict[:outer_shape_bem][:twist][:grid]
     twist_vals = sec_Dict[:outer_shape_bem][:twist][:values]
-    twist_d = FLOWMath.akima(twist_grid,twist_vals,span) .* 180/pi
+    twist_d = safeakima(twist_grid,twist_vals,span) .* 180/pi
 
     chord_grid = sec_Dict[:outer_shape_bem][:chord][:grid]
     chord_vals = sec_Dict[:outer_shape_bem][:chord][:values]
-    chord = FLOWMath.akima(chord_grid,chord_vals,span)
+    chord = safeakima(chord_grid,chord_vals,span)
 
     pitch_axis_grid = sec_Dict[:outer_shape_bem][:pitch_axis][:grid]
     pitch_axis_vals = sec_Dict[:outer_shape_bem][:pitch_axis][:values]
-    pitch_axis = FLOWMath.akima(pitch_axis_grid,pitch_axis_vals,span)
+    pitch_axis = safeakima(pitch_axis_grid,pitch_axis_vals,span)
 
     xoffset = pitch_axis
     aerocenter = pitch_axis #TODO: this was originally used for the automated flutter analysis within the original OWENS code, which is not currently implemented
@@ -145,11 +145,11 @@ function readNuMadGeomCSV(NuMad_geom_file::OrderedCollections.OrderedDict{Symbol
 
             start_nd_arc_grid = layer_Dict[:start_nd_arc][:grid]
             start_nd_arc_vals = layer_Dict[:start_nd_arc][:values]
-            start_nd_arc = FLOWMath.akima(start_nd_arc_grid,start_nd_arc_vals,span)
+            start_nd_arc = safeakima(start_nd_arc_grid,start_nd_arc_vals,span)
 
             end_nd_arc_grid = layer_Dict[:end_nd_arc][:grid]
             end_nd_arc_vals = layer_Dict[:end_nd_arc][:values]
-            end_nd_arc = FLOWMath.akima(end_nd_arc_grid,end_nd_arc_vals,span)
+            end_nd_arc = safeakima(end_nd_arc_grid,end_nd_arc_vals,span)
         else 
             # println("web rotation and offset")
 
@@ -209,16 +209,16 @@ function readNuMadGeomCSV(NuMad_geom_file::OrderedCollections.OrderedDict{Symbol
         # get the number of plies
         n_plies_grid = layer_Dict[:n_plies][:grid]
         n_plies_vals = layer_Dict[:n_plies][:values]
-        stack_layers[:,istack] = FLOWMath.akima(n_plies_grid,n_plies_vals,span) #note that since we cut off layers in the stack sequences below, extrapolated values are not used here
+        stack_layers[:,istack] = safeakima(n_plies_grid,n_plies_vals,span) #note that since we cut off layers in the stack sequences below, extrapolated values are not used here
 
         if !(contains(layer_Dict[:name],"web")) 
             start_nd_arc_grid = layer_Dict[:start_nd_arc][:grid]
             start_nd_arc_vals = layer_Dict[:start_nd_arc][:values]
-            start_nd_arc = FLOWMath.akima(start_nd_arc_grid,start_nd_arc_vals,span)
+            start_nd_arc = safeakima(start_nd_arc_grid,start_nd_arc_vals,span)
 
             end_nd_arc_grid = layer_Dict[:end_nd_arc][:grid]
             end_nd_arc_vals = layer_Dict[:end_nd_arc][:values]
-            end_nd_arc = FLOWMath.akima(end_nd_arc_grid,end_nd_arc_vals,span)
+            end_nd_arc = safeakima(end_nd_arc_grid,end_nd_arc_vals,span)
 
             for ispan = 1:length(span)
                 for iseg=1:length(common_segments)
@@ -322,7 +322,7 @@ end
     #                 println("# midpoint has grid and vals")
     #                 midpoint_nd_arc_grid = layer_Dict[:midpoint_nd_arc][:grid]
     #                 midpoint_nd_arc_vals = layer_Dict[:midpoint_nd_arc][:values]
-    #                 # midpoint_nd_arc = FLOWMath.akima(midpoint_nd_arc_grid,midpoint_nd_arc_vals,span)
+    #                 # midpoint_nd_arc = safeakima(midpoint_nd_arc_grid,midpoint_nd_arc_vals,span)
     #             catch # midpoint is LE or TE
     #                 println("# midpoint is LE or TE")
     #                 if layer_Dict[:midpoint_nd_arc][:fixed] == "LE"
@@ -337,7 +337,7 @@ end
 
     #             width_grid = layer_Dict[:width][:grid]
     #             width_vals = layer_Dict[:width][:values]
-    #             width = FLOWMath.akima(width_grid,width_vals,midpoint_nd_arc_grid)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
+    #             width = safeakima(width_grid,width_vals,midpoint_nd_arc_grid)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
     #             start_nd_arc = midpoint_nd_arc_vals - width./2
     #             end_nd_arc = midpoint_nd_arc_vals + width./2
     #             tryfailed = false
@@ -354,7 +354,7 @@ end
     #                     println("#start had grid and vals")
     #                     start_nd_arc_grid = layer_Dict[:start_nd_arc][:grid]
     #                     start_nd_arc_vals = layer_Dict[:start_nd_arc][:values]
-    #                     # start_nd_arc = FLOWMath.akima(start_nd_arc_grid,start_nd_arc_vals,span)
+    #                     # start_nd_arc = safeakima(start_nd_arc_grid,start_nd_arc_vals,span)
     #                 catch # start is LE or TE
     #                     println("# start is LE or TE")
     #                     if layer_Dict[:start_nd_arc][:fixed] == "LE"
@@ -374,7 +374,7 @@ end
     #                         @error "width grid must be the same span as th start grid, and if the start_nd_arc is fixed, that is 0 to 1"
     #                     end
     #                     width_vals = layer_Dict[:width][:values]
-    #                     width = FLOWMath.akima(width_grid,width_vals,start_nd_arc_grid)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio) to be converted to the arc definition where 0.5 is leading edge, so span is approx 2x chord.
+    #                     width = safeakima(width_grid,width_vals,start_nd_arc_grid)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio) to be converted to the arc definition where 0.5 is leading edge, so span is approx 2x chord.
     #                     end_nd_arc_vals = start_nd_arc_vals+width
 
     #                 catch #otherwise, end must be specified
@@ -382,7 +382,7 @@ end
     #                     try #end has grid and vals
     #                         # end_nd_arc_grid = layer_Dict[:end_nd_arc][:grid]
     #                         end_nd_arc_vals = layer_Dict[:end_nd_arc][:values]
-    #                         # end_nd_arc = FLOWMath.akima(end_nd_arc_grid,end_nd_arc_vals,span)
+    #                         # end_nd_arc = safeakima(end_nd_arc_grid,end_nd_arc_vals,span)
     #                     catch # end is LE or TE
     #                         println("# end is LE or TE")
     #                         if layer_Dict[:end_nd_arc][:fixed] == "LE"
@@ -408,7 +408,7 @@ end
     #                 try #end has grid and vals
     #                     end_nd_arc_grid = layer_Dict[:end_nd_arc][:grid]
     #                     end_nd_arc_vals = layer_Dict[:end_nd_arc][:values]
-    #                     # end_nd_arc = FLOWMath.akima(end_nd_arc_grid,end_nd_arc_vals,span)
+    #                     # end_nd_arc = safeakima(end_nd_arc_grid,end_nd_arc_vals,span)
     #                 catch # end is LE or TE
     #                     println("# end is LE or TE")
     #                     if layer_Dict[:end_nd_arc][:fixed] == "LE"
@@ -426,7 +426,7 @@ end
     #                     @error "width grid must be the same span as th end grid, and if the end_nd_arc is fixed, that is 0 to 1"
     #                 end
     #                 width_vals = layer_Dict[:width][:values]
-    #                 width = FLOWMath.akima(width_grid,width_vals,end_nd_arc_grid)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
+    #                 width = safeakima(width_grid,width_vals,end_nd_arc_grid)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
     #                 start_nd_arc_vals = end_nd_arc_vals-width
 
     #                 tryfailed = false
@@ -443,7 +443,7 @@ end
     #                 @error "offset_y_pa_grid grid must span 0 to 1, set thickness to 0 for the layer for the span positions it is not present"
     #             end
     #             offset_y_pa_vals = layer_Dict[:offset_y_pa][:values]
-    #             offset_y_pa = FLOWMath.akima(offset_y_pa_grid,offset_y_pa_vals,span) # in meters, need to convert to side
+    #             offset_y_pa = safeakima(offset_y_pa_grid,offset_y_pa_vals,span) # in meters, need to convert to side
 
     #             side = layer_Dict[:side]
 
@@ -455,7 +455,7 @@ end
 
     #             width_grid = layer_Dict[:width][:grid]
     #             width_vals = layer_Dict[:width][:values]
-    #             width = FLOWMath.akima(width_grid,width_vals,span)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
+    #             width = safeakima(width_grid,width_vals,span)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
     #             start_nd_arc_vals = midpoint_nd_arc .- width./2
     #             end_nd_arc_vals = midpoint_nd_arc .+ width./2
 
@@ -478,15 +478,15 @@ end
 
     # rotation_grid = layer_Dict[:rotation][:grid]
     # rotation_vals = layer_Dict[:rotation][:values]
-    # rotation = FLOWMath.akima(rotation_grid,rotation_vals,span) # in meters, need to convert to side
+    # rotation = safeakima(rotation_grid,rotation_vals,span) # in meters, need to convert to side
 
     # offset_y_pa_grid = layer_Dict[:offset_y_pa][:grid]
     # offset_y_pa_vals = layer_Dict[:offset_y_pa][:values]
-    # offset_y_pa = FLOWMath.akima(offset_y_pa_grid,offset_y_pa_vals,span) # in meters, need to convert to side
+    # offset_y_pa = safeakima(offset_y_pa_grid,offset_y_pa_vals,span) # in meters, need to convert to side
 
     # width_grid = layer_Dict[:width][:grid]
     # width_vals = layer_Dict[:width][:values]
-    # width = FLOWMath.akima(width_grid,width_vals,span)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
+    # width = safeakima(width_grid,width_vals,span)./(chord*2) # unifying around start and end arc points requires width in meters (as defined by windio to be converted to the arc definition where)
     # start_nd_arc = midpoint_nd_arc - width./2
     # end_nd_arc = midpoint_nd_arc + width./2
 
@@ -700,7 +700,7 @@ function saveOWENSfiles(filename,mymesh,myort,myjoint,myEl,pBC,numadIn_bld)
 
     # Save Blade
     # Used
-    chord = FLOWMath.akima(numadIn_bld.span./maximum(numadIn_bld.span),numadIn_bld.chord,mymesh.structuralSpanLocNorm[1,:])
+    chord = safeakima(numadIn_bld.span./maximum(numadIn_bld.span),numadIn_bld.chord,mymesh.structuralSpanLocNorm[1,:])
 
     bldArray = zeros(length(mymesh.structuralSpanLocNorm),16)
     row = 1
diff --git a/src/meshing_utilities.jl b/src/meshing_utilities.jl
index 55d011f3..3134f00c 100644
--- a/src/meshing_utilities.jl
+++ b/src/meshing_utilities.jl
@@ -37,13 +37,13 @@ function mesh_beam(;L1 = 31.5, #first section of beam length
 
         # First Section
         mesh_x1 = collect(LinRange(zeroOffset,zeroOffset+L1,Nelem1+1))
-        mesh_y1 = FLOWMath.akima(x_shape,y_shape,mesh_x1)
-        mesh_z1 = FLOWMath.akima(x_shape,z_shape,mesh_x1)
+        mesh_y1 = safeakima(x_shape,y_shape,mesh_x1)
+        mesh_z1 = safeakima(x_shape,z_shape,mesh_x1)
 
         # Angled Section
         mesh_x2 = zeroOffset+L1.+collect(LinRange(0.0,L2,Nelem2+1))
-        mesh_y2 = FLOWMath.akima(x_shape,y_shape,mesh_x2)
-        mesh_z2 = FLOWMath.akima(x_shape,z_shape,mesh_x2)
+        mesh_y2 = safeakima(x_shape,y_shape,mesh_x2)
+        mesh_z2 = safeakima(x_shape,z_shape,mesh_x2)
     end
 
     # intra-beam connectivity
@@ -163,13 +163,13 @@ function mesh_beam_centered(;L1 = 6.0, #first section of beam length
 
         # First Section
         mesh_x1 = collect(LinRange(zeroOffset,zeroOffset+L1,Nelem1+1))
-        mesh_y1 = FLOWMath.akima(x_shape,y_shape,mesh_x1)
-        mesh_z1 = FLOWMath.akima(x_shape,z_shape,mesh_x1)
+        mesh_y1 = safeakima(x_shape,y_shape,mesh_x1)
+        mesh_z1 = safeakima(x_shape,z_shape,mesh_x1)
 
         # Angled Section
         mesh_x2 = zeroOffset+L1.+collect(LinRange(0.0,L2,Nelem2+1))
-        mesh_y2 = FLOWMath.akima(x_shape,y_shape,mesh_x2)
-        mesh_z2 = FLOWMath.akima(x_shape,z_shape,mesh_x2)
+        mesh_y2 = safeakima(x_shape,y_shape,mesh_x2)
+        mesh_z2 = safeakima(x_shape,z_shape,mesh_x2)
     end
 
     # Push a shaft onto the mesh
@@ -401,13 +401,13 @@ function create_mesh_struts(;Htwr_base = 15.0,
         # Ensure the blade shape conforms to the turbine height and radius specs
         bshapex = R .* bshapex./maximum(bshapex)
         bshapez = Hbld .* bshapez./maximum(bshapez)
-        bld_Y = FLOWMath.akima(bshapez,bshapex,bld_Z)
+        bld_Y = safeakima(bshapez,bshapex,bld_Z)
     end
 
     if bshapey == zeros(nbelem+1)
         bld_X = zero(bld_Y)
     else
-        bld_X = FLOWMath.akima(bshapez,bshapey,bld_Z)
+        bld_X = safeakima(bshapez,bshapey,bld_Z)
     end
 
     # AeroDyn Compatability
@@ -760,7 +760,7 @@ function create_arcus_mesh(;
         bld_Y = R.*(1.0.-4.0.*(bld_Z/Hbld.-.5).^2)
     else
         # Ensure the blade shape conforms to the turbine height and radius specs
-        bld_Y = FLOWMath.akima(bshapez,bshapex,bld_Z)
+        bld_Y = safeakima(bshapez,bshapex,bld_Z)
     end
     bld_X = zero(bld_Y)
 
@@ -1612,30 +1612,30 @@ function getSectPropsFromOWENSPreComp(usedUnitSpan,numadIn,precompoutput;GX=fals
     # Now create the splines and sample them at the used span
     origUnitSpan = numadIn.span./numadIn.span[end]
     usedUnitSpan = usedUnitSpan./maximum(usedUnitSpan)
-    ei_flap_used = FLOWMath.akima(origUnitSpan,ei_flap,usedUnitSpan)
-    ei_lag_used = FLOWMath.akima(origUnitSpan,ei_lag,usedUnitSpan)
-    gj_used = FLOWMath.akima(origUnitSpan,gj,usedUnitSpan)
-    ea_used = FLOWMath.akima(origUnitSpan,ea,usedUnitSpan)
-    s_fl_used = FLOWMath.akima(origUnitSpan,s_fl,usedUnitSpan)
-    s_af_used = FLOWMath.akima(origUnitSpan,s_af,usedUnitSpan)
-    s_al_used = FLOWMath.akima(origUnitSpan,s_al,usedUnitSpan)
-    s_ft_used = FLOWMath.akima(origUnitSpan,s_ft,usedUnitSpan)
-    s_lt_used = FLOWMath.akima(origUnitSpan,s_lt,usedUnitSpan)
-    s_at_used = FLOWMath.akima(origUnitSpan,s_at,usedUnitSpan)
-    x_sc_used = FLOWMath.akima(origUnitSpan,x_sc,usedUnitSpan)
-    y_sc_used = FLOWMath.akima(origUnitSpan,y_sc,usedUnitSpan)
-    x_tc_used = FLOWMath.akima(origUnitSpan,x_tc,usedUnitSpan)
-    y_tc_used = FLOWMath.akima(origUnitSpan,y_tc,usedUnitSpan)
-    mass_used = FLOWMath.akima(origUnitSpan,mass,usedUnitSpan)
-    flap_iner_used = FLOWMath.akima(origUnitSpan,flap_iner,usedUnitSpan)
-    lag_iner_used = FLOWMath.akima(origUnitSpan,lag_iner,usedUnitSpan)
-    tw_iner_d_used = FLOWMath.akima(origUnitSpan,tw_iner_d,usedUnitSpan)
-    x_cm_used = FLOWMath.akima(origUnitSpan,x_cm,usedUnitSpan).*0.0
-    y_cm_used = FLOWMath.akima(origUnitSpan,y_cm,usedUnitSpan).*0.0
-
-    ac_used = FLOWMath.akima(origUnitSpan,numadIn.aerocenter,usedUnitSpan)
-    twist_d_used = FLOWMath.akima(origUnitSpan,numadIn.twist_d,usedUnitSpan)
-    chord_used = FLOWMath.akima(origUnitSpan,numadIn.chord,usedUnitSpan)
+    ei_flap_used = safeakima(origUnitSpan,ei_flap,usedUnitSpan)
+    ei_lag_used = safeakima(origUnitSpan,ei_lag,usedUnitSpan)
+    gj_used = safeakima(origUnitSpan,gj,usedUnitSpan)
+    ea_used = safeakima(origUnitSpan,ea,usedUnitSpan)
+    s_fl_used = safeakima(origUnitSpan,s_fl,usedUnitSpan)
+    s_af_used = safeakima(origUnitSpan,s_af,usedUnitSpan)
+    s_al_used = safeakima(origUnitSpan,s_al,usedUnitSpan)
+    s_ft_used = safeakima(origUnitSpan,s_ft,usedUnitSpan)
+    s_lt_used = safeakima(origUnitSpan,s_lt,usedUnitSpan)
+    s_at_used = safeakima(origUnitSpan,s_at,usedUnitSpan)
+    x_sc_used = safeakima(origUnitSpan,x_sc,usedUnitSpan)
+    y_sc_used = safeakima(origUnitSpan,y_sc,usedUnitSpan)
+    x_tc_used = safeakima(origUnitSpan,x_tc,usedUnitSpan)
+    y_tc_used = safeakima(origUnitSpan,y_tc,usedUnitSpan)
+    mass_used = safeakima(origUnitSpan,mass,usedUnitSpan)
+    flap_iner_used = safeakima(origUnitSpan,flap_iner,usedUnitSpan)
+    lag_iner_used = safeakima(origUnitSpan,lag_iner,usedUnitSpan)
+    tw_iner_d_used = safeakima(origUnitSpan,tw_iner_d,usedUnitSpan)
+    x_cm_used = safeakima(origUnitSpan,x_cm,usedUnitSpan).*0.0
+    y_cm_used = safeakima(origUnitSpan,y_cm,usedUnitSpan).*0.0
+
+    ac_used = safeakima(origUnitSpan,numadIn.aerocenter,usedUnitSpan)
+    twist_d_used = safeakima(origUnitSpan,numadIn.twist_d,usedUnitSpan)
+    chord_used = safeakima(origUnitSpan,numadIn.chord,usedUnitSpan)
 
     if GX
 
@@ -1730,8 +1730,8 @@ function getSectPropsFromOWENSPreComp(usedUnitSpan,numadIn,precompoutput;GX=fals
             myxpts_top = LinRange(xaf_top[1],xaf_top[end],nspl)
             myxpts_bot = LinRange(xaf_bot[1],xaf_bot[end],nspl)
 
-            myypts_top = FLOWMath.akima(xaf_top,yaf_top,myxpts_top)
-            myypts_bot = FLOWMath.akima(reverse(xaf_bot),reverse(yaf_bot),reverse(myxpts_bot))
+            myypts_top = safeakima(xaf_top,yaf_top,myxpts_top)
+            myypts_bot = safeakima(reverse(xaf_bot),reverse(yaf_bot),reverse(myxpts_bot))
 
             myxafpc[ipci,:] = [myxpts_top;myxpts_bot;myxpts_top[1]]
             myyafpc[ipci,:] = [myypts_top;reverse(myypts_bot);myypts_top[1]]
@@ -1748,11 +1748,11 @@ function getSectPropsFromOWENSPreComp(usedUnitSpan,numadIn,precompoutput;GX=fals
         myxaf = zeros(length(usedUnitSpan)-1,nspl*2+1)
         myyaf = zeros(length(usedUnitSpan)-1,nspl*2+1)
         myzaf = LinRange(0,1,length(usedUnitSpan)-1)
-        chord = FLOWMath.akima(myzafpc,mychord,myzaf)
+        chord = safeakima(myzafpc,mychord,myzaf)
 
         for iline = 1:nspl*2+1
-            myxaf[:,iline] = FLOWMath.akima(myzafpc,myxafpc[:,iline],myzaf)
-            myyaf[:,iline] = FLOWMath.akima(myzafpc,myyafpc[:,iline],myzaf)
+            myxaf[:,iline] = safeakima(myzafpc,myxafpc[:,iline],myzaf)
+            myyaf[:,iline] = safeakima(myzafpc,myyafpc[:,iline],myzaf)
         end
     else
         myxaf = nothing
@@ -1874,7 +1874,7 @@ function create_hawt_mesh(;
     if maximum(bshapez) != 0.0 #TODO more robust
         bshapez = tip_precone .* bshapez./maximum(bshapez)
     end
-    bld_Z = FLOWMath.akima(bshapex,bshapez,bld_Y)
+    bld_Z = safeakima(bshapex,bshapez,bld_Y)
 
     bld_X = zero(bld_Z)
 
@@ -2129,22 +2129,22 @@ function create_hawt_biwing_mesh(;
 
     # Root
     bld_root_Y = collect(LinRange(0.0,R_root,nbelem_root+1))
-    bld_root_Z = FLOWMath.akima(bshapex_root,bshapez_root,bld_root_Y)    
+    bld_root_Z = safeakima(bshapex_root,bshapez_root,bld_root_Y)    
     bld_root_X = zero(bld_root_Y)
 
     # Biwing upper
     bld_biwing_U_Y = collect(LinRange(R_root,R_biwing,nbelem_biwing+1))
-    bld_biwing_U_Z = FLOWMath.akima(bshapex_biwing_U,bshapez_biwing_U,bld_biwing_U_Y)
+    bld_biwing_U_Z = safeakima(bshapex_biwing_U,bshapez_biwing_U,bld_biwing_U_Y)
     bld_biwing_U_X = zero(bld_biwing_U_Y)
 
     # Biwing lower
     bld_biwing_L_Y = collect(LinRange(R_root,R_biwing,nbelem_biwing+1))
-    bld_biwing_L_Z = FLOWMath.akima(bshapex_biwing_L,bshapez_biwing_L,bld_biwing_L_Y)
+    bld_biwing_L_Z = safeakima(bshapex_biwing_L,bshapez_biwing_L,bld_biwing_L_Y)
     bld_biwing_L_X = zero(bld_biwing_L_Y)
 
     # Tip
     bld_tip_Y = collect(LinRange(R_biwing,R_tip,nbelem_tip+1))
-    bld_tip_Z = FLOWMath.akima(bshapex_tip,bshapez_tip,bld_tip_Y)
+    bld_tip_Z = safeakima(bshapex_tip,bshapez_tip,bld_tip_Y)
     bld_tip_X = zero(bld_tip_Y)
 
     # Iterate around the hub
diff --git a/src/utilities.jl b/src/utilities.jl
index 99f27014..b514b05d 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -761,3 +761,11 @@ end
 #     end
 #     return plystress
 # end #stress_calc
+
+function safeakima(x,y,xpt)
+    if minimum(xpt)<(minimum(x)-abs(minimum(x))*0.1) || maximum(xpt)>(maximum(x)+abs(maximum(x))*0.1)
+        msg="Extrapolating on akima spline results in undefined solutions minimum(xpt)<minimum(x) $(minimum(xpt))<$(minimum(x)) or maximum(xpt)<maximum(x) $(maximum(xpt))>$(maximum(x))"
+        throw(OverflowError(msg))
+    end
+    return FLOWMath.akima(x,y,xpt)
+end
\ No newline at end of file
diff --git a/test/Fig4_5_campbell2.jl b/test/Fig4_5_campbell2.jl
index 37b3a2a9..c064078c 100644
--- a/test/Fig4_5_campbell2.jl
+++ b/test/Fig4_5_campbell2.jl
@@ -68,9 +68,9 @@ SNL34m_5_3_Torque = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Data/5.3_T
 
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = OWENS.safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = OWENS.safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 offsetTime = 20.0 # seconds
 tocp = [0.0;new_t.+offsetTime; 1e6]
@@ -86,7 +86,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = OWENS.safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
diff --git a/test/check_mesh.jl b/test/check_mesh.jl
index cf05081d..9c24b600 100644
--- a/test/check_mesh.jl
+++ b/test/check_mesh.jl
@@ -180,34 +180,34 @@ newspan2 = (newspan[1:end-1]+newspan[2:end])/2
 newspan2 = newspan2.-newspan2[1]
 
 tw_aero_new = [sectionPropsArray_bld[ii].twist[1] for ii = 1:length(sectionPropsArray_bld)]
-tw_aero_spl = FLOWMath.akima(newspan2,tw_aero_new,span_loc)
+tw_aero_spl = OWENS.safeakima(newspan2,tw_aero_new,span_loc)
 
 ei_flap_new = [sectionPropsArray_bld[ii].EIyy[1] for ii = 1:length(sectionPropsArray_bld)]
-ei_flap_spl = FLOWMath.akima(newspan2,ei_flap_new,span_loc)
+ei_flap_spl = OWENS.safeakima(newspan2,ei_flap_new,span_loc)
 
 ei_lag_new = [sectionPropsArray_bld[ii].EIzz[1] for ii = 1:length(sectionPropsArray_bld)]
-ei_lag_spl = FLOWMath.akima(newspan2,ei_lag_new,span_loc)
+ei_lag_spl = OWENS.safeakima(newspan2,ei_lag_new,span_loc)
 
 gj_new = [sectionPropsArray_bld[ii].GJ[1] for ii = 1:length(sectionPropsArray_bld)]
-gj_spl = FLOWMath.akima(newspan2,gj_new,span_loc)
+gj_spl = OWENS.safeakima(newspan2,gj_new,span_loc)
 
 ea_new = [sectionPropsArray_bld[ii].EA[1] for ii = 1:length(sectionPropsArray_bld)]
-ea_spl = FLOWMath.akima(newspan2,ea_new,span_loc)
+ea_spl = OWENS.safeakima(newspan2,ea_new,span_loc)
 
 mass_new = [sectionPropsArray_bld[ii].rhoA[1] for ii = 1:length(sectionPropsArray_bld)]
-mass_spl = FLOWMath.akima(newspan2,mass_new,span_loc)
+mass_spl = OWENS.safeakima(newspan2,mass_new,span_loc)
 
 flap_iner_new = [sectionPropsArray_bld[ii].rhoIyy[1] for ii = 1:length(sectionPropsArray_bld)]
-flap_iner_spl = FLOWMath.akima(newspan2,flap_iner_new,span_loc)
+flap_iner_spl = OWENS.safeakima(newspan2,flap_iner_new,span_loc)
 
 lag_iner_new = [sectionPropsArray_bld[ii].rhoIzz[1] for ii = 1:length(sectionPropsArray_bld)]
-lag_iner_spl = FLOWMath.akima(newspan2,lag_iner_new,span_loc)
+lag_iner_spl = OWENS.safeakima(newspan2,lag_iner_new,span_loc)
 
 x_cm_new = [sectionPropsArray_bld[ii].zcm[1] for ii = 1:length(sectionPropsArray_bld)]
-x_cm_spl = FLOWMath.akima(newspan2,x_cm_new,span_loc)
+x_cm_spl = OWENS.safeakima(newspan2,x_cm_new,span_loc)
 
 y_cm_new = [sectionPropsArray_bld[ii].ycm[1] for ii = 1:length(sectionPropsArray_bld)]
-y_cm_spl = FLOWMath.akima(newspan2,y_cm_new,span_loc)
+y_cm_spl = OWENS.safeakima(newspan2,y_cm_new,span_loc)
 
 
 for ii = 1:length(sectionPropsArray_bld)