# Simplified ikeda Fortran code to Python


```python
%matplotlib inline
%load_ext autoreload
%autoreload 2
```


```python
from jupyterthemes import jtplot
jtplot.style(theme='chesterish', context='notebook', ticks=True, grid=False)
```


```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from rolldecayestimators.simplified_ikeda import calculate_roll_damping
from rolldecay.sensitivity import variate_ship, plot_variation, calculate, calculate_variation, _plot_result
import sympy as sp
from rolldecayestimators import symbols, equations
from rolldecayestimators.substitute_dynamic_symbols import lambdify
import rolldecayestimators.ikeda_simple  ## Peter Piehl implementation
```


```python
def plot(OG_div_d = -0.2):
    L_div_B = 6.0
    B_div_d = 4.0
    CB = 0.65
    CMID = 0.98
    PHI = 10
    bBK_div_B = 0.025
    lBK_div_Lpp = 0.2
    
    
    LPP = 300
    Beam = LPP/L_div_B
    DRAFT = Beam/B_div_d
    
    lBK = LPP*lBK_div_Lpp
    bBK =Beam*bBK_div_B
    omega_hat = np.linspace(0.25,1,10)
    OMEGA = omega_hat/(np.sqrt(Beam / (2 * 9.81)))
    
    OG = DRAFT*OG_div_d
    
    B44HAT, BFHAT, BWHAT, BEHAT, BBKHAT, BLHAT = calculate_roll_damping(LPP,Beam,CB,CMID,OG,PHI,lBK,bBK,OMEGA,DRAFT)
    
    fig,ax = plt.subplots()
    fig.set_size_inches(8,5)
    label = """
    Lpp:%0.1f
    Beam:%0.1f
    Draught:%0.1f
    lBK:%0.1f
    bBK:%0.1f
    OG:%0.1f
    """ % (LPP,Beam,DRAFT,lBK,bBK,OG)
    ax.plot(omega_hat,B44HAT, label=label)
    ax.set_title('OG/d:%0.1f' % OG_div_d)
    ax.set_xlabel('$\hat{\omega}$')
    ax.set_ylabel('$\hat{B_{44}}$')
    ylim = ax.get_ylim()
    ax.set_ylim(0,ylim[1])
    ax.set_xlim(0.2,1)
    ax.legend()
    ax.grid(True)
```




```python
plot(OG_div_d=-0.2)
```




```python
plot(OG_div_d=-1.5)
```

### Compare with validation


```python
df_verification_raw = pd.read_csv('ikeda_verification.csv', sep=';', index_col=0)
df_verification_raw.fillna(0,inplace=True)
df_verification_raw
```


```python
df_verification = df_verification_raw.drop(columns='unit')
df_verification = df_verification.transpose()
df_verification
```


```python
LPP = df_verification['LPP']
Breadth = df_verification['Breadth']
Draft = df_verification['Draft']
LBKL = df_verification['LBK/LPP']
BBKB = df_verification['BBK/Breadth']
Tw = df_verification['Tw']
OGd = df_verification['OG/d']
Displacement = df_verification['Displacement']
GM = df_verification['GM']

```


```python
df = pd.DataFrame(index= df_verification.index)
df['lpp'] = LPP
df['beam'] = Breadth
df['DRAFT'] = Draft
df['phi_max'] = 1.0  #[deg]
df['BKL'] = LBKL*LPP
df['BKB'] = BBKB*Breadth
df['omega0'] =  2*np.pi/Tw
df['OG'] = OGd*Draft
rho=1000
df['CB'] = Displacement/rho/(LPP*Breadth*Draft)
df['A0'] = 0.95  # Assumption!
df['kg'] = -df['OG'] + Draft
df['GM'] = GM
df['g'] = 9.81
```


```python
ship=df.iloc[0]
```


```python
ship
```


```python
N=40
changes = np.linspace(1,20,N)
plot_variation(ship=ship, key='phi_max', changes=changes)
```

The relationship between B44 and N coefficient (Bertin) is as follows


```python
N = sp.symbols('N')
omega_hat = sp.symbols('omega_hat')
N_equation = sp.Eq(N, symbols.B_44_hat*sp.pi*symbols.beam*omega_hat/(symbols.GM*symbols.phi_a))
N_equation
```


```python
omega_hat_equation = sp.Eq(omega_hat, symbols.omega*sp.sqrt(symbols.beam/(2*symbols.g)))
omega_hat_equation
```


```python
eqs = [omega_hat_equation,
       N_equation,
      ]
N_equation_omega = sp.Eq(N,sp.solve(eqs, omega_hat, N)[N])
N_equation_omega
```


```python
N_lambda = lambdify(sp.solve(N_equation_omega,N)[0])
```


```python
key='phi_max'
df_variation = variate_ship(ship=ship, key=key, changes=changes)
```


```python
result = calculate_variation(df=df_variation)
result = pd.concat([result, df_variation], axis=1)
```


```python
result.head()
```


```python
result['N'] = N_lambda(B_44_hat=result['B44HAT'], GM=result['GM'], beam=result['beam'], g=result['g'],
                     omega=result['omega0'], phi_a=result['phi_max'])
```


```python
result.plot(y=['N'])
```




```python
for ship_name, ship in df.iterrows():
    df_variation = variate_ship(ship=ship, key=key, changes=changes)
    result = calculate_variation(df=df_variation)
    result = pd.concat([result, df_variation], axis=1)
    
    result['N'] = N_lambda(B_44_hat=result['B44HAT'], GM=result['GM'], beam=result['beam'], g=result['g'],
                     omega=result['omega0'], phi_a=result['phi_max'])
    
    fig,ax=plt.subplots()
    result.plot(y=['N'], ax=ax)
    ax.set_title(ship_name)
    
    
    
```

## Froude scaling ok?


```python
df2 = df.copy()
scale=50.0

df2['lpp']*=scale
df2['beam']*=scale
df2['DRAFT']*=scale
#df2['phi_max'] = 1.0  #[deg]
df2['BKL']*=scale
df2['BKB']*=scale
df2['omega0']*=1/np.sqrt(scale)
df2['OG']*=scale
#df2['CB'] = Displacement/rho/(LPP*Breadth*Draft)
#df2['A0'] = 0.95  # Assumption!
df2['kg']*=scale
df2['GM']*=scale
```


```python
for ship_name, ship in df2.iterrows():
    df_variation = variate_ship(ship=ship, key=key, changes=changes)
    result = calculate_variation(df=df_variation)
    result = pd.concat([result, df_variation], axis=1)
    
    result['N'] = N_lambda(B_44_hat=result['B44HAT'], GM=result['GM'], beam=result['beam'], g=result['g'],
                     omega=result['omega0'], phi_a=result['phi_max'])
    
    fig,ax=plt.subplots()
    result.plot(y=['N'], ax=ax)
    ax.set_title(ship_name)
```


```python
model = df.iloc[0]
model
```


```python
ship = df2.iloc[0]
ship
```


```python
N=40
changes = np.linspace(1,20,N)

df_variation = variate_ship(ship=model, key=key, changes=changes)
result_model = calculate_variation(df=df_variation)

df_variation = variate_ship(ship=ship, key=key, changes=changes)
result_ship = calculate_variation(df=df_variation)

result = pd.merge(left=result_model, right=result_ship, how='left', left_index=True, right_index=True, 
                 suffixes=('_m','_s'))

fig,ax=plt.subplots()
result.iloc[:,0:5].plot(ax=ax)
result.iloc[:,5:].plot(ax=ax, style='--')
ax.set_title('Model scale and full scale')

```

...it seems that the friction cannot be neglected in model scale

## Compare with Peter Piehl


```python
OG_div_d = 0
L_div_B = 6.0
B_div_d = 4.0
CB = 0.65
CMID = 0.98
PHI = 10
bBK_div_B = 0.025
lBK_div_Lpp = 0.2


LPP = 300
Beam = LPP/L_div_B
DRAFT = Beam/B_div_d

lBK = LPP*lBK_div_Lpp
bBK =Beam*bBK_div_B
omega_hat = np.linspace(0.25,1,10)
OMEGA = omega_hat/(np.sqrt(Beam / (2 * 9.81)))

OG = DRAFT*OG_div_d

B44HAT, BFHAT, BWHAT, BEHAT, BBKHAT, BLHAT = calculate_roll_damping(LPP,Beam,CB,CMID,OG,PHI,lBK,bBK,OMEGA,DRAFT)
```


```python

```


```python
volume = 164000  # [m3]
#DRAFT = 14.0,  # [m]

data = {
    'lpp'  : 350.0,  # [m]
    #'beam' : 51.0,  # [m]
    'CB'   : 0.655,  # [1]
    'A0' : 0.987,  # [1]
    'CoG'  : 0,  #[m]
    'PHI'  : 20.0,  # [deg]
    'Tw'   : 19.0, # [s]
    'BKL'  : 5*14.85,  # [m]
    'BKB'  : 0.6,  # [m]
}
data['omega0']=2*np.pi/data['Tw']
BdivT=np.linspace(2.55,4.48,20)
df_input = pd.DataFrame(data=data,index=BdivT)
df_input.index.name=r'B/T'
#df_input['DRAFT']=df_input['beam']/df_input.index
df_input['DRAFT']=np.sqrt(volume/(df_input['lpp']*df_input.index*df_input['CB']))
df_input['beam']=volume/(df_input['lpp']*df_input['DRAFT']*df_input['CB'])

df_input['kg']=df_input['DRAFT']+df_input['CoG']
df_input['phi_max']=df_input['PHI']  # Deg!
df_input
```


```python
df_result = calculate_variation(df_input)
df_result
```


```python
fig,ax=plt.subplots()
df_result.plot(ax=ax)
ax.grid(True)
ax.set_xlim(2.5,4.5)
ax.set_ylim(0,0.018)
```




```python
ikeda = rolldecayestimators.ikeda_simple.Ikeda()
```


```python
para = {\
        'LPP':370.0,\
        'BRTH':51.0,\
        'DRAFT':14.0,\
        'CB':0.655,\
        'CMID':0.987,\
        'OG':0,\
        'PHI':20.0,\
        'TW':19.0,\
        'BKCOMP':True,\
        'lBK':5*14.85,\
        'bBK':0.6,\
        }
```


```python
ikeda.setPara(para)
ikeda.ikedaMethod()
BHAT = ikeda.getBHAT()
```


```python
renames = {
    'lpp':'LPP',
    'beam':'BRTH',
    'A0':'CMID',
    'CoG':'OG',
    'TW':'Tw',
    'BKL':'lBK',
    'BKB':'bBK',
}
for old_key,new_key in renames.items():
    df_input[new_key]=df_input[old_key]

data = df_input.iloc[0]

```


```python
ikeda.setPara(data)
ikeda.ikedaMethod()
BHAT = ikeda.getBHAT()
```


```python
def calculate_variation_piehl(row):
    
    row=row.copy()
    row['BKCOMP'] = (row['BKL']>0)
    
    ikeda=rolldecayestimators.ikeda_simple.Ikeda()
    ikeda.setPara(row)
    ikeda.ikedaMethod()
    BHAT = ikeda.getBHAT()
    result = pd.Series(name=row.name)
    result['B44HAT'] = ikeda.B44HAT
    result['BFHAT'] = ikeda.BFHAT
    result['BWHAT'] = ikeda.BWHAT
    result['BEHAT'] = ikeda.BEHAT
    result['BBKHAT'] = ikeda.BBKHAT
    return result
```


```python
df_result2 = df_input.apply(func=calculate_variation_piehl, axis=1)
```


```python
fig,ax=plt.subplots()
df_result.plot(ax=ax, alpha=0.6)
df_result2.plot(ax=ax, style='--')
ax.grid(True)
ax.set_xlim(2.5,4.5)
ax.set_ylim(0,0.018)
```


```python

```

8333e4c45c3da47409eaabb79e79a1787ec71596
ipynb
Jupyter Notebook
notebooks/rolldecay/04_simplified_ikeda/06.1_maa_simplified_ikeda_fortran.ipynb
martinlarsalbert/rolldecay
bee335c27e6519bd32ed9488751bc13e7f6fa9b3
MIT
notebooks/rolldecay/04_simplified_ikeda/06.1_maa_simplified_ikeda_fortran.ipynb
martinlarsalbert/rolldecay
bee335c27e6519bd32ed9488751bc13e7f6fa9b3
MIT
2021-02-02T23:07:50.000Z
2022-01-13T03:27:41.000Z
notebooks/rolldecay/04_simplified_ikeda/06.1_maa_simplified_ikeda_fortran.ipynb
martinlarsalbert/rolldecay
bee335c27e6519bd32ed9488751bc13e7f6fa9b3
MIT
Qwen/Qwen-72B
1. YES
2. YES