Jovian
Sign In
3 years ago
#! /usr/bin/env python3

import numpy as np
#from sincfit import *
from scipy.optimize import root
import matplotlib.pyplot as plt
import matplotlib as mpl

def yalpha(alpha):
    if alpha==0:
        return 1
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        toint=1-1/sqas*np.arctan(sqas)
        integ=np.trapz(toint,svals)
        return np.real(np.sqrt(1+alpha)-integ/np.pi)

def yaprime(alpha):
    dalpha = 1e-6
    if np.abs(alpha) < dalpha:
        return 1/6
    
    return np.real(0.5*(yalpha(alpha+dalpha)-yalpha(alpha-dalpha))/dalpha)

def yadiff(alpha):
    if alpha==0:
        return 1/6
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        ass=alpha/(1+svals**2)
        # sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        sqas=np.lib.scimath.sqrt(ass)
        
        toint=1/sqas*np.arctan(sqas)-1/(1+ass)
        integ=np.trapz(toint,svals)
        return np.real(1/(2*np.sqrt(1+alpha))-integ/(2*np.pi*alpha))

def flux0(B,fu):
    #phi_0 = (chi_0 * f.u.)/(mu_b ( mu_0)) * B for Bexternal
    chi0 = 0.274 #external
    flux0 = (fu/(4*np.pi*9.274*0.1))*chi0*B
    return flux0

def kappaout(kappain,uguess,T,TMF,T0,Q,J,flux0):
    kB = 0.0861733 # Boltzmann constant in meV/K
    kperp = kappain[0]
    kpara = kappain[1]
    alpha = (kperp-kpara)/(kperp+kpara)
    gamma = uguess*kB*Q/J
    y = yalpha(alpha)
    yprime = yadiff(alpha)
    k23 = 1/((kperp+kpara)**(3/2))
    kperpline = (T-TMF)/T0 + 1/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(3/2*1/np.sqrt(kperp+kpara)*y+(4*kpara-2*kperp)*k23*yprime )-kperp
    kparaline = (T-TMF)/T0 + 3/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(2*1/np.sqrt(kperp+kpara)*y + (2*kpara-6*kperp)*k23*yprime )-kpara
    ksquared = np.array((kperpline,kparaline))
    return ksquared

T0 = 1.13 # K
TMF = 30.6 # K
Q = 0.039 # A-1
J = 2.8 # meV/A
fu = 24.018 # A3
kperp=0.01 # A-1
kpara=0.02 # A-1

kappain = np.array(((kperp/Q)**2,(kpara/Q)**2))

uguess = 229
#uguess = 127

B = 0.104

# for B in (0.104,0.240,0.340,0.350):
phi0 = flux0(B,fu)
for T in np.arange(27.0,32.0,0.1):
    #print(B,T)
    if B==0.0:
        kpara=0.010 # A-1
        kperp=0.010 # A-1
        
    huh = root(kappaout,kappain,args=(uguess,T,TMF,T0,Q,J,phi0))
    # print(huh)

    if(huh.success):
        kout=np.sqrt((huh.x))*Q
        print(B, T, *kout)
c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:55: RuntimeWarning: invalid value encountered in double_scalars c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:56: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:57: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:39: RuntimeWarning: invalid value encountered in sqrt
0.104 29.00000000000003 0.0058998727469731765 0.017851499993875117 0.104 29.10000000000003 0.006714095941307619 0.01858899428166738 0.104 29.20000000000003 0.007652430587795579 0.01935300204552879 0.104 29.300000000000033 0.008718437764471697 0.020151657107193847 0.104 29.400000000000034 0.00991198183834818 0.0209952871236127 0.104 29.500000000000036 0.011229026429167366 0.021895208426311326 0.104 29.600000000000037 0.012661574150001985 0.022861973584683226 0.104 29.70000000000004 0.01419776135210306 0.02390342139341465 0.104 29.80000000000004 0.015822236243081943 0.025023029597085772 0.104 29.90000000000004 0.017516974397525864 0.026219060404547554 0.104 30.000000000000043 0.01926254800083737 0.027484762416842404 0.104 30.100000000000044 0.021039641142570398 0.02880952183690363 0.104 30.200000000000045 0.022830462843858323 0.030180537345192323 0.104 30.300000000000047 0.024619754338055022 0.031584500483421825 0.104 30.40000000000005 0.02639526626579444 0.03300890341175285 0.104 30.50000000000005 0.028147764372642606 0.03444283178428824 0.104 30.60000000000005 0.02987072051918108 0.03587729072454394 0.104 30.700000000000053 0.031559853822473574 0.03730520055413953 0.104 30.800000000000054 0.03321264482678747 0.0387212050163696 0.104 30.900000000000055 0.03482789367067327 0.04012140102518611 0.104 31.000000000000057 0.036405352202353525 0.041503058122005596 0.104 31.10000000000006 0.03794543465695998 0.042864363062477755 0.104 31.20000000000006 0.039448998901079696 0.04420420354824746 0.104 31.30000000000006 0.04091718572121294 0.04552199309789779 0.104 31.400000000000063 0.04235130338882501 0.04681753329658288 0.104 31.500000000000064 0.04375274635044801 0.04809090755705519 0.104 31.600000000000065 0.045122939045013806 0.04934240030915814 0.104 31.700000000000067 0.04646329791255937 0.050572436138280996 0.104 31.800000000000068 0.047775206413481305 0.05178153432712241 0.104 31.90000000000007 0.04905999924563856 0.052970275120073906 
#! /usr/bin/env python3

import numpy as np
#from sincfit import *
from scipy.optimize import root
import matplotlib.pyplot as plt
import matplotlib as mpl

def yalpha(alpha):
    if alpha==0:
        return 1
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        toint=1-1/sqas*np.arctan(sqas)
        integ=np.trapz(toint,svals)
        return np.real(np.sqrt(1+alpha)-integ/np.pi)

def yaprime(alpha):
    dalpha = 1e-6
    if alpha < dalpha:
        return 1/6
    
    return np.real(0.5*(yalpha(alpha+dalpha)-yalpha(alpha-dalpha))/dalpha)

def yadiff(alpha):
    if alpha==0:
        return 1/6
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        ass=alpha/(1+svals**2)
        # sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        sqas=np.lib.scimath.sqrt(ass)
        
        toint=1/sqas*np.arctan(sqas)-1/(1+ass)
        integ=np.trapz(toint,svals)
        return np.real(1/(2*np.sqrt(1+alpha))-integ/(2*np.pi*alpha))

def flux0(B,fu):
    #phi_0 = (chi_0 * f.u.)/(mu_b ( mu_0)) * B for Bexternal
    chi0 = 0.274 #external
    flux0 = (fu/(4*np.pi*9.274*0.1))*chi0*B
    return flux0

def kappaout(kappain,uguess,T,TMF,T0,Q,J,flux0):
    kB = 0.0861733 # Boltzmann constant in meV/K
    kperp = kappain[0]
    kpara = kappain[1]
    alpha = (kperp-kpara)/(kperp+kpara)
    gamma = uguess*kB*Q/J
    y = yalpha(alpha)
    yprime = yadiff(alpha)
    k23 = 1/((kperp+kpara)**(3/2))
    kperpline = (T-TMF)/T0 + 1/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(3/2*1/np.sqrt(kperp+kpara)*y+(4*kpara-2*kperp)*k23*yprime )-kperp
    kparaline = (T-TMF)/T0 + 3/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(2*1/np.sqrt(kperp+kpara)*y + (2*kpara-6*kperp)*k23*yprime )-kpara
    ksquared = np.array((kperpline,kparaline))
    return ksquared

T0 = 1.13 # K
TMF = 30.6 # K
Q = 0.039 # A-1
J = 2.8 # meV/A
fu = 24.018 # A3
kperp=0.01 # A-1
kpara=0.02 # A-1

kappain = np.array(((kperp/Q)**2,(kpara/Q)**2))

uguess = 85 #5964
#uguess = 127

B = 0.240

# for B in (0.104,0.240,0.340,0.350):
phi0 = flux0(B,fu)
for T in np.arange(27.0,32.0,0.1):
    #print(B,T)
    if B==0.0:
        kpara=0.10 # A-1
        kperp=0.10 # A-1
        
    huh = root(kappaout,kappain,args=(uguess,T,TMF,T0,Q,J,phi0))
    # print(huh)
    if(huh.success):
        kout=np.sqrt((huh.x))*Q
        print(B, T, *kout)
c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:55: RuntimeWarning: invalid value encountered in double_scalars c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:56: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:57: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:39: RuntimeWarning: invalid value encountered in sqrt
0.24 29.90000000000004 0.004778436060905078 0.0276172798196949 0.24 30.000000000000043 0.00860769236673203 0.02888546497583155 0.24 30.100000000000044 0.012121068081599086 0.030301525889325 0.24 30.200000000000045 0.015328018919181344 0.031819472387130045 0.24 30.300000000000047 0.01826632157848297 0.033391534055159294 0.24 30.40000000000005 0.020975372738500007 0.03498209987152442 0.24 30.50000000000005 0.023489945549917336 0.03656786899199603 0.24 30.60000000000005 0.025839119904827085 0.03813467615892608 0.24 30.700000000000053 0.028046732253682666 0.03967442666854008 0.24 30.800000000000054 0.03013221842705006 0.04118286750104785 0.24 30.900000000000055 0.032111458786495746 0.04265811717310517 0.24 31.000000000000057 0.03399750364874809 0.044099733368102315 0.24 31.10000000000006 0.035801158690873275 0.04550813060436817 0.24 31.20000000000006 0.03753144497742358 0.04688421836889386 0.24 31.30000000000006 0.03919595681805712 0.048229176610198904 0.24 31.400000000000063 0.040801139696751504 0.04954431674536336 0.24 31.500000000000064 0.04235250677202574 0.050830996096976964 0.24 31.600000000000065 0.04385480845739237 0.052090565878667526 0.24 31.700000000000067 0.045312166162224804 0.05332434033885748 0.24 31.800000000000068 0.046728178558994554 0.05453357928250187 0.24 31.90000000000007 0.04810600667315999 0.055719479046833435 
#! /usr/bin/env python3

import numpy as np
#from sincfit import *
from scipy.optimize import root
import matplotlib.pyplot as plt
import matplotlib as mpl

def yalpha(alpha):
    if alpha==0:
        return 1
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        toint=1-1/sqas*np.arctan(sqas)
        integ=np.trapz(toint,svals)
        return np.real(np.sqrt(1+alpha)-integ/np.pi)

def yaprime(alpha):
    dalpha = 1e-6
    if np.abs(alpha) < dalpha:
        return 1/6
    
    return np.real(0.5*(yalpha(alpha+dalpha)-yalpha(alpha-dalpha))/dalpha)

def yadiff(alpha):
    if alpha==0:
        return 1/6
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        ass=alpha/(1+svals**2)
        # sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        sqas=np.lib.scimath.sqrt(ass)
        
        toint=1/sqas*np.arctan(sqas)-1/(1+ass)
        integ=np.trapz(toint,svals)
        return np.real(1/(2*np.sqrt(1+alpha))-integ/(2*np.pi*alpha))

def flux0(B,fu):
    #phi_0 = (chi_0 * f.u.)/(mu_b ( mu_0)) * B for Bexternal
    chi0 = 0.274 #external
    flux0 = (fu/(4*np.pi*9.274*0.1))*chi0*B
    return flux0

def kappaout(kappain,uguess,T,TMF,T0,Q,J,flux0):
    kB = 0.0861733 # Boltzmann constant in meV/K
    kperp = kappain[0]
    kpara = kappain[1]
    alpha = (kperp-kpara)/(kperp+kpara)
    gamma = uguess*kB*Q/J
    y = yalpha(alpha)
    yprime = yadiff(alpha)
    k23 = 1/((kperp+kpara)**(3/2))
    kperpline = (T-TMF)/T0 + 1/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(3/2*1/np.sqrt(kperp+kpara)*y+(4*kpara-2*kperp)*k23*yprime )-kperp
    kparaline = (T-TMF)/T0 + 3/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(2*1/np.sqrt(kperp+kpara)*y + (2*kpara-6*kperp)*k23*yprime )-kpara
    ksquared = np.array((kperpline,kparaline))
    return ksquared

T0 = 1.13 # K
TMF = 30.6 # K
Q = 0.039 # A-1
J = 2.8 # meV/A
fu = 24.018 # A3
kperp=0.02 # A-1
kpara=0.05 # A-1

kappain = np.array(((kperp/Q)**2,(kpara/Q)**2))

uguess = 13419 #5964
#uguess = 127

B = 0.350

# for B in (0.104,0.240,0.340,0.350):
phi0 = flux0(B,fu)
for T in np.arange(27.0,32.0,0.1):
    #print(B,T)
    if B==0.0:
        kpara=0.01 # A-1
        kperp=0.01 # A-1
        
    huh = root(kappaout,kappain,args=(uguess,T,TMF,T0,Q,J,phi0))
    # print(huh)
    if(huh.success):
        kout=np.sqrt((huh.x))*Q
        print(B, T, *kout)
c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in sqrt c:\mcstas-2.6\miniconda3\lib\site-packages\ipykernel_launcher.py:39: RuntimeWarning: invalid value encountered in sqrt 
#! /usr/bin/env python3

import numpy as np
#from sincfit import *
from scipy.optimize import root
import matplotlib.pyplot as plt
import matplotlib as mpl

def yalpha(alpha):
    if alpha==0:
        return 1
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        toint=1-1/sqas*np.arctan(sqas)
        integ=np.trapz(toint,svals)
        return np.real(np.sqrt(1+alpha)-integ/np.pi)

def yaprime(alpha):
    dalpha = 1e-6
    if alpha < dalpha:
        return 1/6
    
    return np.real(0.5*(yalpha(alpha+dalpha)-yalpha(alpha-dalpha))/dalpha)

def yadiff(alpha):
    if alpha==0:
        return 1/6
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        ass=alpha/(1+svals**2)
        # sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        sqas=np.lib.scimath.sqrt(ass)
        
        toint=1/sqas*np.arctan(sqas)-1/(1+ass)
        integ=np.trapz(toint,svals)
        return np.real(1/(2*np.sqrt(1+alpha))-integ/(2*np.pi*alpha))

def flux0(B,fu):
    #phi_0 = (chi_0 * f.u.)/(mu_b ( mu_0)) * B for Bexternal
    chi0 = 0.274 #external
    flux0 = (fu/(4*np.pi*9.274*0.1))*chi0*B
    return flux0

def kappaout(kappain,uguess,T,TMF,T0,Q,J,flux0):
    kB = 0.0861733 # Boltzmann constant in meV/K
    kperp = kappain[0]
    kpara = kappain[1]
    alpha = (kperp-kpara)/(kperp+kpara)
    gamma = uguess*kB*Q/J
    y = yalpha(alpha)
    yprime = yadiff(alpha)
    k23 = 1/((kperp+kpara)**(2/3))
    kperpline = (T-TMF)/T0 + 1/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(3/2*1/np.sqrt(kperp+kpara)*y+(4*kpara-2*kperp)*k23*yprime )-kperp
    kparaline = (T-TMF)/T0 + 3/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(2*1/np.sqrt(kperp+kpara)*y + (2*kpara-6*kperp)*k23*yprime )-kpara
    ksquared = np.array((kperpline,kparaline))
    return ksquared

T0 = 1.13 # K
TMF = 30.6 # K
Q = 0.039 # A-1
J = 2.8 # meV/A
fu = 24.018 # A3
kperp=0.020 # A-1
kpara=0.050 # A-1

kappain = np.array(((kperp/Q)**2,(kpara/Q)**2))

uguess = 5964
#uguess = 127

B = 0.384

# for B in (0.104,0.240,0.340,0.350):
phi0 = flux0(B,fu)
for T in np.arange(27.0,32.0,0.1):
    #print(B,T)
    if B==0.0:
        kpara=0.01 # A-1
        kperp=0.01 # A-1
        
    huh = root(kappaout,kappain,args=(uguess,T,TMF,T0,Q,J,phi0))
    # print(huh)
    if(huh.success):
        kout=np.sqrt((huh.x))*Q
        print(B, T, *kout)
0.384 27.0 0.30515592524783164 0.45131851591699973 0.384 27.1 0.3054913428885933 0.45140555554247397 0.384 27.200000000000003 0.3058261598584991 0.45149235707336566 0.384 27.300000000000004 0.30616037875313 0.4515789211965764 0.384 27.400000000000006 0.30649400214885797 0.4516652485948612 0.384 27.500000000000007 0.3068270326030434 0.4517513399468646 0.384 27.60000000000001 0.3071594726542292 0.4518371959271571 0.384 27.70000000000001 0.30749132482233343 0.4519228172062698 0.384 27.80000000000001 0.30782259160883774 0.4520082044507305 0.384 27.900000000000013 0.3081532754969764 0.4520933583230972 0.384 28.000000000000014 0.30848337895191796 0.4521782794819942 0.384 28.100000000000016 0.30881290442094866 0.45226296858214443 0.384 28.200000000000017 0.30914185433365404 0.4523474262744036 0.384 28.30000000000002 0.3094702311020943 0.45243165320579315 0.384 28.40000000000002 0.30979803712097964 0.4525156500195334 0.384 28.50000000000002 0.31012527476784413 0.45259941735507603 0.384 28.600000000000023 0.3104519464032168 0.452682955848135 0.384 28.700000000000024 0.31077805437078904 0.4527662661307196 0.384 28.800000000000026 0.3111036009975816 0.4528493488311652 0.384 28.900000000000027 0.3114285885941101 0.45293220457416394 0.384 29.00000000000003 0.3117530194545454 0.4530148339807959 0.384 29.10000000000003 0.3120768958568755 0.453097237668559 0.384 29.20000000000003 0.31240022006306284 0.4531794162513991 0.384 29.300000000000033 0.31272299431920164 0.45326137033973923 0.384 29.400000000000034 0.3130452208556709 0.4533431005405096 0.384 29.500000000000036 0.31336690188728755 0.4534246074571759 0.384 29.600000000000037 0.31368803961345754 0.45350589168976824 0.384 29.70000000000004 0.31400863621832364 0.45358695383490955 0.384 29.80000000000004 0.3143286938709127 0.45366779448584316 0.384 29.90000000000004 0.3146482147252806 0.45374841423246093 0.384 30.000000000000043 0.31496720092065555 0.4538288136613307 0.384 30.100000000000044 0.31528565458157964 0.453908993355723 0.384 30.200000000000045 0.3156035778173047 0.4539889538958782 0.384 30.300000000000047 0.3159209727248349 0.454068695858098 0.384 30.40000000000005 0.31623784138480704 0.4541482198161386 0.384 30.50000000000005 0.3165541858643989 0.45422752634034796 0.384 30.60000000000005 0.31687000821678446 0.454306615997908 0.384 30.700000000000053 0.31718531048126575 0.4543854893528599 0.384 30.800000000000054 0.31750009468340185 0.45446414696612963 0.384 30.900000000000055 0.3178143628351379 0.454542589395553 0.384 31.000000000000057 0.3181281169349326 0.45462081719589936 0.384 31.10000000000006 0.31844135896788106 0.4546988309188974 0.384 31.20000000000006 0.31875409090584206 0.4547766311132578 0.384 31.30000000000006 0.3190663147075582 0.4548542183246976 0.384 31.400000000000063 0.31937803231877776 0.4549315930959638 0.384 31.500000000000064 0.3196892456723733 0.45500875596685725 0.384 31.600000000000065 0.3199999566884615 0.45508570747425453 0.384 31.700000000000067 0.3203101672745181 0.4551624481521317 0.384 31.800000000000068 0.32061987932549413 0.45523897853158696 0.384 31.90000000000007 0.3209290947239299 0.4553152991408623 
#! /usr/bin/env python3

import numpy as np
#from sincfit import *
from scipy.optimize import root
import matplotlib.pyplot as plt
import matplotlib as mpl

def yalpha(alpha):
    if alpha==0:
        return 1
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        toint=1-1/sqas*np.arctan(sqas)
        integ=np.trapz(toint,svals)
        return np.real(np.sqrt(1+alpha)-integ/np.pi)

def yaprime(alpha):
    dalpha = 1e-6
    if np.abs(alpha) < dalpha:
        return 1/6
    
    return np.real(0.5*(yalpha(alpha+dalpha)-yalpha(alpha-dalpha))/dalpha)

def yadiff(alpha):
    if alpha==0:
        return 1/6
    
    else:
        svals=np.linspace(-10000,10000,10000000,endpoint=True)
        ass=alpha/(1+svals**2)
        # sqas=np.lib.scimath.sqrt(alpha/(1+svals**2))
        sqas=np.lib.scimath.sqrt(ass)
        
        toint=1/sqas*np.arctan(sqas)-1/(1+ass)
        integ=np.trapz(toint,svals)
        return np.real(1/(2*np.sqrt(1+alpha))-integ/(2*np.pi*alpha))

def flux0(B,fu):
    #phi_0 = (chi_0 * f.u.)/(mu_b ( mu_0)) * B for Bexternal
    chi0 = 0.274 #external
    flux0 = (fu/(4*np.pi*9.274*0.1))*chi0*B
    return flux0

def kappaout(kappain,uguess,T,TMF,T0,Q,J,flux0):
    kB = 0.0861733 # Boltzmann constant in meV/K
    kperp = kappain[0]
    kpara = kappain[1]
    alpha = (kperp-kpara)/(kperp+kpara)
    gamma = uguess*kB*Q/J
    y = yalpha(alpha)
    yprime = yadiff(alpha)
    k23 = 1/((kperp+kpara)**(3/2)) # not 2/3
    kperpline = (T-TMF)/T0 + 1/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(3/2*1/np.sqrt(kperp+kpara)*y+(4*kpara-2*kperp)*k23*yprime ) #-kperp
    kparaline = (T-TMF)/T0 + 3/6*uguess*flux0**2 + 1/6*1/(np.sqrt(2)*np.pi)*gamma*T*(2*1/np.sqrt(kperp+kpara)*y + (2*kpara-6*kperp)*k23*yprime ) #-kpara
    ksquared = np.array((kperpline,kparaline))
    return ksquared

T0 = 1.13 # K
TMF = 30.6 # K
Q = 0.039 # A-1
J = 2.8 # meV/A
fu = 24.018 # A3
kperp=0.01319 # A-1
kpara=0.02607 # A-1
T = 30
B = 0.240
#30	0.013188540636044	0.0260761677385981
#kperp=0.01 #A-1
#kpara=0.01 #A-1
#T = 30.0
#B = 0.0

kappain = np.array(((kperp/Q)**2,(kpara/Q)**2))

u = 0.25 # meV/A3

ucalc = u/(J*Q*Q)

uguess = 85.0 #229.0 #314.0# 5964
#uguess = 96.1

phi0 = flux0(B,fu)

#kappasquared = kappaout(kappain,uguess,T,TMF,T0,Q,J,phi0)
#print(kappasquared)

alpha = (kperp**2-kpara**2)/(kperp**2+kpara**2)

huh = kappaout(kappain,uguess,T,TMF,T0,Q,J,phi0)
#print(huh)
kout=np.sqrt((huh))*Q

print('B = {:.3f} T'.format(B))
print('T = {:.1f} K'.format(T))
print('~u = {:.2f} calculated'.format(ucalc))
print('~u = {:.2f} guessed'.format(uguess))
print('alpha = {:.4f}'.format(alpha))
print('y(a) = {:.4f}'.format(yalpha(alpha)))
print('y\'(a) = {:.4f}'.format(yaprime(alpha)))
print('input:')
print('k_perp = {:.5f} A-1, k_para = {:.5f} A-1'.format(kperp,kpara))
print('output:')
print('k_perp = {:.5f} A-1, k_para = {:.5f} A-1'.format(*kout))
B = 0.240 T T = 30.0 K ~u = 58.70 calculated ~u = 85.00 guessed alpha = -0.5924 y(a) = 0.8898 y'(a) = 0.2122 input: k_perp = 0.01319 A-1, k_para = 0.02607 A-1 output: k_perp = 0.00605 A-1, k_para = 0.02851 A-1