# This program calculates J'/J, carrying errors through correctly
# follows methodology from http://stacks.iop.org/1367-2630/10/i=8/a=083025

from uncertainties import ufloat

rho_over_J = 0.183
b          = 2.43
mu_B       = ufloat((9.27400915e-24,0.000023e-26)) # http://physics.nist.gov/cgi-bin/cuu/Value?mub
k_B        = ufloat((1.3806504e-23,0.0000024e-23)) # http://physics.nist.gov/cgi-bin/cuu/Value?k

input_filename  = 'beta-vs-JJ-in.tsv'
output_filename = 'beta-vs-JJ.tsv'

import csc
from uncertainties import unumpy as unp
#import numpy as np

pi = 3.14159265358979323

title,input,error = csc.read(input_filename)

#output = csc.begin({'title':'beta vs J\'/J'},['beta','beta_err','JJ','JJ_err','#name'],output_filename)

for compound in input:
    g    = ufloat((float(compound['g']), float(compound['g_err']) ))
    Bc   = ufloat((float(compound['Bc']),float(compound['Bc_err'])))
    TN   = ufloat((float(compound['TN']),float(compound['TN_err'])))
    
    # work out the J'/J value
    J             = (mu_B*g*Bc)/(4*k_B)
    Jprime_over_J = unp.exp(b-4*pi*rho_over_J*J/TN)
    
    # pass beta through
    beta = compound['beta']
    beta_err = compound['beta_err']

    # modify the name
    name = '#'+compound['name'] # add a # at the start to comment it out for GLE
    
    print name, J
    
    # [()]s are required because unp spits out 0d arrays at the moment...bah
    #csc.append([beta,beta_err,Jprime_over_J[()].nominal_value,Jprime_over_J[()].std_dev(),name],output)

#csc.close(output)