######################################################################################
######################################################################################
#           Early-stage ship powering prediction (see DDS 051-1, MAY 1982)           #
#                     Code written by Jay Borthen, April 2013                        #
######################################################################################
######################################################################################

# A:  Section area
# A_fins:  Area of fin centerplane times the number of fins
# ATTC:  American Towing Tank Conference
# A_v:  Ship above-water transverse area
# A_x:  Section area at Station of Maximum Area
# A_0:  Section area at Station 0
# A_20:  Section area at Station 20
# BWL:  Maximum beam of waterline at a particular draft
# B_x:  Beam at Design Waterline and at Station of Maximum Area
# B_20:  Beam at Design Waterline and at Station 20
# C_A:  Correlation-allowance coefficient.  This is a function of bottom paint and
#       the condition of the bottom (e.g. roughness, fouling, etc.)
# C_AA:  Air drag coefficient based on ship frontal area
# C_D:  Drag coefficient
# C_D_Ap:  Appendage drag coefficient.  Because the formula herein for P_EAp is not non-dimensional, the values
#         of C_DAp which is used to compute P_EAp are applicable for English-unit computations only.
# C_DBTD:  Drag coefficient for bow thruster duct openings
# C_DFin:  Drag coefficient for stabilizer fins
# C_F:  Frictional resistance coefficient.  Tabulated values given in "Coefficients
#       for International Towing Tank Conference 1957 Model Ship Correlation Line"
#       authored by Hadler, J.B., April 1958
# C_P:  Longitudinal prismatic coefficient
# C_R:  Residuary resistance coefficient
# CRPP:  Controllable-reversible-pitch propeller
# C_RTSS:  Residuary resistance coefficient for TSS hull forms
# C_S:  Wetted surface coefficient.  C_S = (S / (Nabla*LWL)^0.5) or alternatively determined from Figure 2 in DDS 051-1.
# C_STSS:  Wetted surface coefficient for TSS hull forms
# C_WP:  Waterplane area coefficient at Design Waterline
# C_StillAir:  Still-air drag coefficient
# C_x:  Section coefficient at station of maximum area (i.e. Midship coefficient)
# C_Nabla:  Volumetric coefficient
# Delta:  Ship displacement
# D_BTD:  Diameter of bow thruster duct openings
# D_P:  Propeller diameter
# DWL:  Design waterline
# EAR:  Propeller Expanded-Area Ratio
# FPP:  Fixed-Pitch propeller
# ITTC:  International Towing Tank Conference
# i_E:  Design waterline entrance half-angle in degrees
# J:  Propeller advance coefficient
# K_T:  Propeller thrust coefficient
# K_l:  Constant required to convert to standard units of power
# L:  Ship length (usually LBP)
# L_PP:  Length between perpendiculars
# LWL:  Length on waterline at Design Waterline
# PD:  Propeller pitch-to-diameter ratio
# P_E:  Effective power
# P_EAA:  Effective power due to still-air drag
# P_EAP:  Effective power due to appendages.
# P_EBH:  Bare hull effective power
# P_EBTD:  Effective power due to bow thruster duct openings
# P_EFin:  Effective power due to stabilizer fins
# P_EMisc:  Effective power due to miscellaneous appendage items, hull openings, etc.
# P_ETot:  Ship total effective power, inclusive of effective power added by still-air drag and by power margin
# PMF:  Power margin factor.  This "shall be applied to effective power, over the
#       entire speed range" (perr DDS 051-1).
# P_S:  Shaft power
# R_AA:  Resistance due to still air
# R_AP:  Resistance of appendages
# R_BH:  Bare hull resistance
# R_F:  Ship frictional resistance
# Rn:  Reynolds number
# R_R:  Residuary resistance
# R_RShip:  Ship residuary resistance for a specific hull form
# R_RTSS:  Residuary resistance of equivalent Taylor Standard Series hull form
# R_RPerTon:  Residuary resistance per ton of displacement, R_RPerTon = R_R / Delta
# R_T:  Total resistance, R_T = R_BH + R_AP + R_AA
# rps:  revolutions per second
# S:  Wetted surface
# S_Ship:  Wetted surface of specific ship
# S_TSS:  Wetted surface of equivalent TSS hull form.  Determined from Figure
#         2 in DDS 051-1 (MAY 1984).  The figure included in DDS 051-1 is from
#         "Systematic Resistance Experiments with Taylor Models Having a Beam-
#         Draft Ratio of 4.50" authored by Hahnel, G. et al., July 1966
# SWSF: Ship wetted surface factor (i.e. S_Ship/S_TSS).  This is a function
#       of general ship type, bow bulb or dome, transom size and skeg size.
#       The bow dome and skeg wetted surface is typically included in SWSF.
# SWSF_alt:  Alternate ship wetted surface factor (e.g. S/((Delta*LWL)^0.5) where S is in ft^2, Delta is in
#            long tons, and LWL is in ft).  This is typically a function of C_p, C_x, Delta/((LWL/100)^3), 
#            B/T, general ship type, bow bulb, and transom size.  Use Figure 3 from DDS 051-1 to assess 
#            effects of C_x and B/T.
# t:  Thrust-deduction fraction, t = (Thrust - R_T) / Thrust
# T_Mean:  Ship mean draft
# TSS:  Taylor Standard Series
# T_x:  Ship draft on design waterline at station of maximum area
# TWL:  Ship draft to a particular waterline
# V_A:  Propeller speed of advance
# V_S:  Ship speed
# WCF:  Worm curve factor.  Function of speed range.
# w:  Taylor wake fraction, w = (V_S - V_A) / V_S
# Delta_App:  Displacement of appendages
# Delta_BareHull:  Bare hull displacement to Design Waterline
# Delta_Total:  Displacement of fully-appended hull to Design Waterline
# Nabla:  Ship submerged volume
# Nabla_BareHull:  Bare hull volume to Design Waterline
# eta_D:  Propulsive coefficient
# eta_O:  Propeller Open-water efficiency
# eta_R:  Propeller relative-rotative efficiency
# lambda:  Ship model linear scale ratio
# nu:  Kinematic viscosity of salt water at specified temperature and 3.5% salinity
# rho:  Sea water mass density, 1.9905 slugs/ft^2 (to convert from Table 4 in DDS 051-1 values to g/cm^3, 
#       multiply by 0.01602 and then by 32.174)
# rho_A:  Mass density of air

# designStage:
#
#   Stage 1:  During Feasibility and Preliminary Design, prior to development of 
#             a preliminary body plan, appendage configuration, etc.
#   Stage 2:  During Preliminary and Contract Design, prior to conduct of self-
#             propelled model tests
#   Stage 3:  During Preliminary and Contract Design, after self-propelled model
#             tests with the stock propeller have been conducted.
#   Stage 4:  During the final stages of Contract Design, after self-propelled
#             model tests with the design propeller have been conducted.

# hullCondition:  For U.S. Navy ships having the Navy vinyl paint system applied
#                 over sandblasted bottom plating, and having been out of drydock
#                 for only a short time, hullCondition=0. For estimating the 
#                 resistance of U.S. Navy ships two years after the initial 
#                 drydocking, hullCondition=1.

# KinematicViscosityTable:  Table of salt water kinematic viscosity values at 
#                           different temperatures. 


######################################################################################
######################################################################################



# Initialize variables
C_A<-0


# =====================
# Ship Characteristics:
# =====================

ShipType<-3  # 1 = Aircraft Carriers, 2 = Destroyer-Type Ships, 3 = Naval Auxiliaries
LWL<-465.88  # in ft
B_x<-62.00  # in ft
T_x<-20.00  # in ft
Delta_BareHull<-8400  # in long tons
Delta_App<-100 # in long tons
Delta_Total<-Delta_BareHull+Delta_App
C_x<-0.825
C_p<-0.610
BT<-B_x/T_x  # Beam to draft ratio
Nabla_BareHull<-294000  # in ft^3
C_Nabla_BareHull<-0.002908
D_P<-17  # in ft
S<-31849 # in ft^2   <-- if "S" is unknown, approximate via S = 1.025*Lpp*(C_B*B+1.7*T) = 1.025*((Nabla_BareHull/T_x)+1.7*Lpp*T_x)
A_v<-4030  # in ft^2



# ================================
# Constraints, Ratios, and Factors:
# ================================
Temp_SW<-59  # in degrees F
rho<-1.9905  # in  (lb)*(sec^2)*(ft^4)
nu<-0.000012817  # in (ft^2)/(sec)

# C_A_1 is calculated via formula from "Prediction of Resistance and Propulsion Power of Ships" (Kristensen, et al, 2012)
C_A_1<-(0.5*log10(Delta_Total*1.01605)-0.1*(log10(Delta_Total*1.01605))^2)/1000  # 1.01605 is conversion factor from long tons to tonnes

C_A<-0.0004  <-- Typical default value
# if (LWL>=190 & LWL<=960){
#   C_A<-(0.008289/(LWL^(1/3)))-0.00064
# }else if (LWL<190){
#   C_A<-0.0008
# }else{
#   C_A<-0.0002
# }
# if (hullCondition==1){
#   C_A<-C_A+0.0007
# }

C_S_TSS<-2.54  # Visually estimated from DDS 051-1 Figure 2
C_S_Ship = (S/((Nabla_BareHull*LWL)^0.5))

(C_S_Ship/C_S_TSS)

(S/((Delta_Total*LWL)^0.5))

PropNum<-2  # shafts assumed to be strut-supported

if (PropNum==2){
  
  C_D_AP<-((0.000004*LWL^2)-(0.0066*LWL)+5.0053)*10^(-5)
  
}else if (PropNum==1){
  
  # TBD
  
}else{
  
  paste("Number of props is outside bounds.")
  
}


if (ShipType==1){
  C_AA<-0.45
}else if (ShipType==2){
  C_AA<-0.70
}else if (ShipType==3){
  C_AA<-0.75  
}else{
  paste("Ship type unknown.")
}

hullCondition<-0
designStage<-1

VL<-seq(0.60,1.3,by=0.1)  # EXAMPLE
V_knots<-VL*(LWL^0.5)  # EXAMPLE
WCF<-c(3.45,3.22,2.20,1.54,1.09,0.92,0.83,0.80)  # EXAMPLE

#VL<-seq(0.50,0.95,by=0.05)  # VL: Speed-to-Length^(1/2) ratio
#V_knots<-VL*(LWL^0.5)
#WCF_AO177<-c(0.89,0.87,0.91,0.92,0.88,0.80,0.73,0.81,0.95,1.03)


Fn<-0.2976*V_knots/sqrt(LWL)  # V_knots in knots, LWL in ft, and g in ft/sec^2.  0.2976 = (1.68781 ft/sec/knot)/sqrt(32.174 ft/sec^2)

if (designStage==0){
  margin<-10
}else if (designStage==1){
  margin<-8
}else if (designStage==2){
  margin<-6
}else if (designStage==3){
  margin<-4
}
PMF<-(1+(margin/100))

# Calculate Reynold's number:
#KinematicViscosityTable<-read.table("SaltWaterKinematicViscosityTable.csv",sep=",",header=TRUE)
#nu<-KinematicViscosityTable[which(((KinematicViscosityTable[,1])==waterTemp)),2]  # in (m^2/s)*10^6
Rn<-(V_knots*1.68781*LWL)/nu  # 1.68781 is the conversion factor from knots to ft/sec

# Determine frictional resistance coefficient:
C_F<-0.075/((log10(Rn)-2)^2)  # Utilizes ITTC Line

# Calculate ship frictional resistance:
#WaterDensityTable<-read.table("WaterDensityTable.csv",sep=",",header=TRUE)
#rho<-WaterDensityTable[which(((WaterDensityTable[,1])==waterTemp)),2]  # in (g/m^3)
#rho<-rho*0.01602*32.174
R_F<-(rho/2)*S*((V_knots*1.68781)^2)*(C_F+C_A)  # 1.68781 is the conversion factor from knots to ft/sec
S_TSS<-C_S_TSS*(Nabla_BareHull*LWL)^0.5
LNabla<-round(LWL/(Nabla_BareHull^(1/3)),1)  # LNabla: Length-to-Volume^(1/3) ratio


# The following C_R equation is from Kristensen, et al.
C_R<-(40733.90490*(Fn^5)-36751.17452*(Fn^4)+13395.20097*(Fn^3)-2415.43170*(Fn^2)+215.01129*(Fn)-7.26)/1000

# At the same value of V/sqrt(LWL), R_RTSS/Delta_TSS = R_RShip/Delta_Ship
# R_RTSS_to_Delta<- rho*S_TSS*(V_knots^2)*C_R/(2*Delta)
R_R_TSS<- rho*S_TSS*(V_knots^2)*C_R/2

R_R_Ship<-R_R_TSS*WCF

R_T<-R_F+R_R

C_D_Ap<-(-1.862*10^-09)*LWL^3+(4.615*10^-6)*LWL^2-(4.306*10^-3)*LWL+3.536  # Equation from MNVDET

EHP<-(C_T*rho*S*V_knots^3)/(2*550)  # 550 ft-lb/sec