function heat_capacity_nb3sn(c_pt,T,B,TC,TCS,TC0) result(heat_capacity) bind(C)
type(c_ptr), value :: c_pt
real(dp), value :: T,B,TC,TCS,TC0
real(dp) :: heat_capacity
real(dp) :: CP,CPN,CPS,AAA,BBB,CCC,DDD,F,CPN2
real(dp) :: AA,BB,CC,DD,a,bbbb,c,d,na,nb,nc,nd,TMAX,TT
type(nb3sn_cfg_t), pointer :: nb3sn_ptr
call c_f_pointer(c_pt,nb3sn_ptr)
AA = 38.2226876_dp ; BB = -848.36422_dp
CC = 1415.13807_dp ; DD = -346.83796_dp
a = 6.804586085_dp ; bbbb = 59.92091818_dp
c = 25.82863336_dp ; d = 8.779183354_dp
na = 1.0_dp ; nb = 2.0_dp
nc = 3.0_dp ; nd = 4.0_dp
TMAX = 400.0_dp
if(T<=TC) then
TT=TC/TC0
CCC=(-0.46306_dp) - (0.067830_dp)*TC
DDD=27.2_dp/(1.0_dp+(0.34_dp*TT))**2
AAA=1500.0_dp*(CCC**2)/(2.0_dp*DDD-1.0_dp)
if (TC<=10.0_dp) then
BBB=(7.5475E-3_dp)*TC**2
else if(TC>10.0_dp .and. TC<=20.0_dp) then
BBB=(-0.3_dp+0.00375_dp*TC**2)/0.09937_dp
endif
CPS=(AAA+BBB)*(T/TC)**3
else
CPS=0.0_dp
endif
if (T<=10.0_dp .and. Tc<=10.0_dp) then
CPN=(7.5475e-3_dp)*T**2
else if (T>10.0_dp .and. T<=20.0_dp .and. Tc>10.0_dp .and. Tc<=20.0_dp) then
CPN=(-0.3_dp + 0.00375_dp*T**2)/(0.09937_dp)
else
TT=T
TT=min(TT,TMAX)
CPN = AA*TT /(a+TT)**na + BB*TT**2/(bbbb+TT)**nb + CC*TT**3/(c+TT)**nc + DD*TT**4/(d+TT)**nd
endif
if (T<=10.0_dp) then
CPN2=(7.5475e-3_dp)*T**2
else if (T>10.0_dp .and. T<=20.0_dp) then
CPN2=(-0.3_dp + 0.00375_dp*T**2)/(0.09937_dp)
else
TT=T
TT=min(TT,TMAX)
CPN2 = AA*TT /(a+TT)**na + BB*TT**2/(bbbb+TT)**nb + CC*TT**3/(c+TT)**nc + DD*TT**4/(d+TT)**nd
endif
if(T<=TCS) then
CP=CPS
else if (T>TCS .and. T<=TC) then
if(TCS<TC) then
F= (T-TCS)/(TC-TCS)
else
F= 1.0_dp
endif
CP= F*CPN + (1.0_dp-F)*CPS
else if (T>TC) then
CP=CPN2
endif
heat_capacity=CP
end function heat_capacity_nb3sn