#! /usr/bin/python3
# Last edited on 2024-07-26 15:36:37 by stolfi

from math import sin, cos, tan, atan2, asin, hypot, sqrt, pi, inf, floor, sqrt
import rn
import sys

def profile(N,nr,ni,nt,th,alpha):
  
  vr = profile_root()
  vi = profile_involute()
  vt = profile_tip()
  
  vp = []
  append_section(vp, vr, +1, 'r0')
  append_section(vp, vi, +1, 'i0')
  append_section(vp, vt, -1, 't0')
  append_section(vp, vt, +1, 't1')
  append_section(vp, vi, -1, 'i1')
  append_section(vp, vi, -1, 'r1')
  return vp
  
def profile_root():
  ctr = 
  rad = 
  t0 = 
  t1 = 
  return circle(ctr,rad,t0,t1,nr)
    
def profile_involute():
  ctr = 
  rad = 
  t0 = 
  t1 = 
  return involute(ctr,rad,t0,t1,nr)

def profile_tip():
  ctr = 
  rad = 
  t0 = 
  t1 = 
  return circle(ctr,rad,t0,t1,nr)
  
def circle(ctr,rad,t0,t1,n):
  vs = []
  for k in range(n+1)
    fk = k/n
    tk = (1-fk)*t0 + fk*t1
    pk = rn.add(ctr, rn.scale(rad, (cos(tk), sin(tk))))
    vs.append(pk)
  return vs
  
def append_section(vp, vs, dir, lab):
  ns = len(vs)-1;
  for k in range(ns+1):
    js = k if dir > 0 else ns-k
    vsj = vs[js]
    vp.append((vsj[0], vsj[1], vsj[2], f"{lab}.{js:04d}"))

def involute(N,n):
  alpha0 = 20*pi/180 # Pressure angle (radians).
  p_c = th           # Circular pitch.
  if N == inf:
    p_a = 0
    d_p = inf
    m = p_c/pi
  else:
    m = p_c/pi      # Modulus. 
    d_p = p_c*N     # Pitch diameter.
    p_a = 2*pi/N    # Angular pitch, radians. 
  alpha = 20 # Pressure angle.
  

  vs = []
  for k in range(n+1)
    d = 
    fk = k/n
    tk = (1-fk)*t0 + fk*t1
    pk = rn.add(ctr, rn.scale(rad, (cos(tk), sin(tk))))
    vs.append(pk)
  return vs
  
 
def inv(a):
    return math.tan(a) - a

def alpha_at_given_delta(d_p,delta,d_b):
    return math.acos(d_b/(d_p+delta))

def thickness_at_given_delta((d_p+delta),d_p,d_b,alpha,m,x):
    S_p = math.pi/2*m + 2*m*x*math.tan(alpha)
    alpha_d = alpha_at_given_d((d_p+delta),d_b)
    S = (d_p+delta)*(S_p/d_p+inv(alpha)-inv(alpha_d))
    return S