from scipy.special import j0
from math import cos, pi

def q_ppp(r,z,P=3000.0E6,E_d=0.9,H=3.8,a=0.005,n=57800.0,E_r=1.0,R_tilde=1.8,
          H_tilde=3.8):
    ''' q_ppp calculates the power density in a reactor given coordinate (r,z)
    
    q_ppp uses the derivation for heat production in fuel elements from Lamarsh
    and Baratta to determine the power density at a given location (r,z) within
    a reactor, given the quantities power, locally deposited energy, 
    recoverable energy, core radius and height including extrapolation, radius
    of the fuel rod, and number of rods.  The keyword arguments are:
    
    Keyword Arguments:
    P -- the power given in Watts (default 3000.0E6 W)
    E_d -- the locally deposited energy per fission (default 0.9)
    E_r -- the recoverable energy per fission in same units as E_d (default 1.0)
    H -- the height of the core in cm (default 380)
    a -- the radius of a rod in cm (default 0.5)
    n -- the number of rods in the core (default 57800)
    R_tilde -- the core radius with extrapolated distance in cm (default 180)
    H_tilde -- the core height with extrapolated distance in cm (default H)
    '''
    q_ppp = 1.16*P*E_d*j0(2.405*r/R_tilde)*cos(pi*z/H_tilde)/(H*a*a*n*E_r);
    units = "W/m^3";
    
    return (q_ppp,units);
    
def T_s(T_m=2200,k_f=3.0,P=3000.0E6,E_d=0.9,H=3.8,a=0.005,n=57800.0,E_r=1.0,
        R_tilde=1.8,H_tilde=3.8):
    ''' T_s finds the surface temperature outside of an ideal fuel rod
    
    T_s uses the resistance heating model to determine the T_s using the 
    assumption that the heat generated by the entire rod is equal to the power
    density times the volume
    
    Keyword Arguments:
    T_m -- the centerline temperature given in degress Celsius (default 2200.0)
    k_f -- the conduction coefficient of the fuel (default 3 W/mK)
    P -- the power given in Watts (default 3000.0E6 W)
    E_d -- the locally deposited energy per fission (default 0.9)
    E_r -- the recoverable energy per fission in same units as E_d (default 1.0)
    H -- the height of the core in cm (default 380)
    a -- the radius of a rod in cm (default 0.5)
    n -- the number of rods in the core (default 57800)
    R_tilde -- the core radius with extrapolated distance in cm (default 180)
    H_tilde -- the core height with extrapolated distance in cm (default H)
    '''
    q_ppp_cl,units = q_ppp(0.0,0.0);
    q = pi*a*a*H*q_ppp_cl;
    R_f = 1/(4*pi*H*k_f);
    T_s = T_m - q*R_f;
    return (T_s,'deg C');
    
def T_m(k_f=3.0,H=3.8):
    return (1/(4*pi*H*k_f))