C*********************************************************************************
C*                                                                               *
C*       First principle equation of state of hydrogen-helium mixture            *
C*       Mixing ratio: 9 helium atoms per 110 hydrogen atoms                     *
C*                                                                               *
C*       Free energy interpolation code written by B. Militzer                   *
C*                                                                               *
C*       Copyright by B. Militzer                       Berkeley, CA 01-29-13    *
C*                                                                               *
C*********************************************************************************

        real*8 FUNCTION F3_58114(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:17)
        real*8 T(1:17)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.002218737897d0 ,
     3  0.003833723896d0 ,
     4  0.005997547323d0 ,
     5  0.008896725291d0 ,
     6  0.01278116095d0 ,
     7  0.01798568464d0 ,
     8  0.02495891555d0 ,
     9  0.03430193178d0 ,
     1  0.04682008205d0 ,
     2  0.06359240505d0 ,
     3  0.08606464039d0 ,
     4  0.1161738484d0 ,
     5  0.1565153805d0 ,
     6  0.2105665936d0 ,
     7  0.2829865905d0 ,
     8  0.3800178256d0 /
        data F /
     1  -0.1084358657d0 ,
     2  -0.1145921372d0 ,
     3  -0.124025543d0 ,
     4  -0.1381498883d0 ,
     5  -0.159005403d0 ,
     6  -0.189550851d0 ,
     7  -0.2343363721d0 ,
     8  -0.3002486779d0 ,
     9  -0.3986073784d0 ,
     1  -0.5414962946d0 ,
     2  -0.7465596719d0 ,
     3  -1.039421767d0 ,
     4  -1.462122235d0 ,
     5  -2.067812907d0 ,
     6  -2.934545472d0 ,
     7  -4.164694056d0 ,
     8  -5.894257843d0 /

        FT1 = -4.858567699d0
        FTN = -18.25028466d0
        FF = splint2(T,F,17,FT1,FTN,TT)
        F3_58114 = FF
        RETURN
        END

        real*8 FUNCTION Fnn3_58114(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:17)
        real*8 T(1:17)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.002218737897d0 ,
     3  0.003833723896d0 ,
     4  0.005997547323d0 ,
     5  0.008896725291d0 ,
     6  0.01278116095d0 ,
     7  0.01798568464d0 ,
     8  0.02495891555d0 ,
     9  0.03430193178d0 ,
     1  0.04682008205d0 ,
     2  0.06359240505d0 ,
     3  0.08606464039d0 ,
     4  0.1161738484d0 ,
     5  0.1565153805d0 ,
     6  0.2105665936d0 ,
     7  0.2829865905d0 ,
     8  0.3800178256d0 /
        data F /
     1  0.1701651336d0 ,
     2  0.345718794d0 ,
     3  0.5751106695d0 ,
     4  0.8722163456d0 ,
     5  1.213129264d0 ,
     6  1.662326535d0 ,
     7  2.309223227d0 ,
     8  3.352591718d0 ,
     9  4.99754489d0 ,
     1  7.483745583d0 ,
     2  10.88502357d0 ,
     3  15.62728918d0 ,
     4  22.9621055d0 ,
     5  34.04788025d0 ,
     6  51.06643245d0 ,
     7  75.59768636d0 ,
     8  109.1352712d0 /

        FT1 = 146.4492454d0
        FTN = 346.4683429d0
        FF = splint2(T,F,17,FT1,FTN,TT)
        Fnn3_58114 = FF
        RETURN
        END

        real*8 FUNCTION F3_05139(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:17)
        real*8 T(1:17)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.002218737897d0 ,
     3  0.003833723896d0 ,
     4  0.005997547323d0 ,
     5  0.008896725291d0 ,
     6  0.01278116095d0 ,
     7  0.01798568464d0 ,
     8  0.02495891555d0 ,
     9  0.03430193178d0 ,
     1  0.04682008205d0 ,
     2  0.06359240505d0 ,
     3  0.08606464039d0 ,
     4  0.1161738484d0 ,
     5  0.1565153805d0 ,
     6  0.2105665936d0 ,
     7  0.2829865905d0 ,
     8  0.3800178256d0 /
        data F /
     1  -0.1078827323d0 ,
     2  -0.1135118781d0 ,
     3  -0.1225489902d0 ,
     4  -0.1358410656d0 ,
     5  -0.1557405533d0 ,
     6  -0.1850728855d0 ,
     7  -0.2281984776d0 ,
     8  -0.2916754519d0 ,
     9  -0.3861526082d0 ,
     1  -0.5229449749d0 ,
     2  -0.7190945528d0 ,
     3  -0.999228991d0 ,
     4  -1.402171487d0 ,
     5  -1.978887924d0 ,
     6  -2.80141874d0 ,
     7  -3.967869347d0 ,
     8  -5.611237801d0 /

        FT1 = -4.520997471d0
        FTN = -17.32003547d0
        FF = splint2(T,F,17,FT1,FTN,TT)
        F3_05139 = FF
        RETURN
        END

        real*8 FUNCTION F2_6(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.1067962321d0 ,
     2  -0.1165752892d0 ,
     3  -0.1360914226d0 ,
     4  -0.1731153515d0 ,
     5  -0.2424362465d0 ,
     6  -0.3728813408d0 ,
     7  -0.6079726765d0 ,
     8  -1.023246473d0 ,
     9  -1.760085218d0 ,
     1  -3.056222446d0 ,
     2  -5.31706341d0 /

        FT1 = -3.632790712d0
        FTN = -16.38244493d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F2_6 = FF
        RETURN
        END

        real*8 FUNCTION F2_3(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.1049587458d0 ,
     2  -0.1138834469d0 ,
     3  -0.1322332308d0 ,
     4  -0.1674251142d0 ,
     5  -0.2346414485d0 ,
     6  -0.3611005424d0 ,
     7  -0.5877012258d0 ,
     8  -0.9871483726d0 ,
     9  -1.692149499d0 ,
     1  -2.928554662d0 ,
     2  -5.088435825d0 /

        FT1 = -3.329599107d0
        FTN = -15.67223476d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F2_3 = FF
        RETURN
        END

        real*8 FUNCTION F2_2(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.1039138015d0 ,
     2  -0.1123069863d0 ,
     3  -0.1302167047d0 ,
     4  -0.164745144d0 ,
     5  -0.231378953d0 ,
     6  -0.3562996828d0 ,
     7  -0.5796787469d0 ,
     8  -0.9733005551d0 ,
     9  -1.667173294d0 ,
     1  -2.882044787d0 ,
     2  -5.005096648d0 /

        FT1 = -3.276086424d0
        FTN = -15.40956101d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F2_2 = FF
        RETURN
        END

        real*8 FUNCTION F2_1(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.102562956d0 ,
     2  -0.1104074116d0 ,
     3  -0.1276891461d0 ,
     4  -0.1615494134d0 ,
     5  -0.2275410576d0 ,
     6  -0.3508338478d0 ,
     7  -0.5708217701d0 ,
     8  -0.9583204864d0 ,
     9  -1.640690126d0 ,
     1  -2.833809603d0 ,
     2  -4.91926237d0 /

        FT1 = -2.987251454d0
        FTN = -15.14831008d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F2_1 = FF
        RETURN
        END

        real*8 FUNCTION F2(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.1006883598d0 ,
     2  -0.1079300637d0 ,
     3  -0.1245336049d0 ,
     4  -0.157722379d0 ,
     5  -0.2231249095d0 ,
     6  -0.3445736759d0 ,
     7  -0.5608718563d0 ,
     8  -0.9418356092d0 ,
     9  -1.612608341d0 ,
     1  -2.783162164d0 ,
     2  -4.829090845d0 /

        FT1 = -2.617563685d0
        FTN = -14.86895136d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F2 = FF
        RETURN
        END

        real*8 FUNCTION F1_86(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.09695054015d0 ,
     2  -0.1031129684d0 ,
     3  -0.1185522382d0 ,
     4  -0.1509336843d0 ,
     5  -0.2153913745d0 ,
     6  -0.3338173046d0 ,
     7  -0.5446721428d0 ,
     8  -0.9158888461d0 ,
     9  -1.568713633d0 ,
     1  -2.707154971d0 ,
     2  -4.697522482d0 /

        FT1 = -1.934443034d0
        FTN = -14.48244664d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F1_86 = FF
        RETURN
        END

        real*8 FUNCTION F1_75(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.09239416576d0 ,
     2  -0.09765745294d0 ,
     3  -0.1122233923d0 ,
     4  -0.1441196277d0 ,
     5  -0.2076459355d0 ,
     6  -0.3232219605d0 ,
     7  -0.5292868827d0 ,
     8  -0.8921765255d0 ,
     9  -1.530563358d0 ,
     1  -2.643450723d0 ,
     2  -4.586718469d0 /

        FT1 = -0.867941558d0
        FTN = -14.15753293d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F1_75 = FF
        RETURN
        END

        real*8 FUNCTION F1_6(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:11)
        real*8 T(1:11)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.003130579022d0 ,
     3  0.006511588575d0 ,
     4  0.01191081187d0 ,
     5  0.02053297391d0 ,
     6  0.03430193178d0 ,
     7  0.05628994364d0 ,
     8  0.09140317949d0 ,
     9  0.1474764316d0 ,
     1  0.2370212863d0 ,
     2  0.3800178256d0 /
        data F /
     1  -0.08274799448d0 ,
     2  -0.0867862302d0 ,
     3  -0.1001299327d0 ,
     4  -0.1317116078d0 ,
     5  -0.1928699981d0 ,
     6  -0.304005328d0 ,
     7  -0.502868062d0 ,
     8  -0.8535173564d0 ,
     9  -1.469753537d0 ,
     1  -2.545793386d0 ,
     2  -4.42732388d0 /

        FT1 = -0.09997203068d0
        FTN = -13.70830285d0
        FF = splint2(T,F,11,FT1,FTN,TT)
        F1_6 = FF
        RETURN
        END

        real*8 FUNCTION F1_5(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  -0.07241583934d0 ,
     2  -0.08241537936d0 ,
     3  -0.1128894112d0 ,
     4  -0.1654885214d0 ,
     5  -0.248253423d0 ,
     6  -0.3758063235d0 ,
     7  -0.569080297d0 ,
     8  -0.8586663057d0 ,
     9  -1.545444533d0 ,
     1  -4.308384924d0 /

        FT1 = -0.6256526348d0
        FTN = -13.35843934d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F1_5 = FF
        RETURN
        END

        real*8 FUNCTION F1_4(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  -0.05774402096d0 ,
     2  -0.06758274672d0 ,
     3  -0.09707104051d0 ,
     4  -0.1473320314d0 ,
     5  -0.2270227804d0 ,
     6  -0.3504760928d0 ,
     7  -0.5380142504d0 ,
     8  -0.8195398323d0 ,
     9  -1.487807229d0 ,
     1  -4.179141642d0 /

        FT1 = -0.2061052303d0
        FTN = -13.01535284d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F1_4 = FF
        RETURN
        END

        real*8 FUNCTION F1_3(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  -0.03713591796d0 ,
     2  -0.04647184987d0 ,
     3  -0.07438471074d0 ,
     4  -0.1220743947d0 ,
     5  -0.1983913539d0 ,
     6  -0.3172312929d0 ,
     7  -0.4985216742d0 ,
     8  -0.7714004106d0 ,
     9  -1.420060956d0 ,
     1  -4.03503121d0 /

        FT1 = -0.07965548847d0
        FTN = -12.64509617d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F1_3 = FF
        RETURN
        END

        real*8 FUNCTION F1_2(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  -0.007673167264d0 ,
     2  -0.01558039772d0 ,
     3  -0.04139312168d0 ,
     4  -0.08614329274d0 ,
     5  -0.1587690591d0 ,
     6  -0.2726031086d0 ,
     7  -0.4470803018d0 ,
     8  -0.7105232901d0 ,
     9  -1.338225196d0 ,
     1  -3.872604344d0 /

        FT1 = 0.06086714913d0
        FTN = -12.2540856d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F1_2 = FF
        RETURN
        END

        real*8 FUNCTION F1_1(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  0.03638246613d0 ,
     2  0.03023946049d0 ,
     3  0.007126550235d0 ,
     4  -0.03433031493d0 ,
     5  -0.1027318801d0 ,
     6  -0.2110341263d0 ,
     7  -0.3780983505d0 ,
     8  -0.6313284016d0 ,
     9  -1.236589714d0 ,
     1  -3.685542112d0 /

        FT1 = 0.1728053813d0
        FTN = -11.84172379d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F1_1 = FF
        RETURN
        END

        real*8 FUNCTION F0_95(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  0.151193021d0 ,
     2  0.1476688709d0 ,
     3  0.1295736889d0 ,
     4  0.09363978205d0 ,
     5  0.03249345795d0 ,
     6  -0.06635798831d0 ,
     7  -0.2205517389d0 ,
     8  -0.4564339428d0 ,
     9  -1.0242408d0 ,
     1  -3.33482534d0 /

        FT1 = -0.5528675441d0
        FTN = -11.18756363d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F0_95 = FF
        RETURN
        END

        real*8 FUNCTION F0_8(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  0.3851557858d0 ,
     2  0.3786041134d0 ,
     3  0.3679905565d0 ,
     4  0.3386328532d0 ,
     5  0.2854737687d0 ,
     6  0.1977439794d0 ,
     7  0.0581740385d0 ,
     8  -0.1577256319d0 ,
     9  -0.6823246644d0 ,
     1  -2.839476354d0 /

        FT1 = -8.394054271d0
        FTN = -10.49077162d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F0_8 = FF
        RETURN
        END

        real*8 FUNCTION F0_7(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  0.6659293127d0 ,
     2  0.6583443873d0 ,
     3  0.6533683214d0 ,
     4  0.629041328d0 ,
     5  0.5828526145d0 ,
     6  0.5046773717d0 ,
     7  0.3783466645d0 ,
     8  0.1789622497d0 ,
     9  -0.3117967737d0 ,
     1  -2.35195104d0 /

        FT1 = -13.09918801d0
        FTN = -9.959013659d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F0_7 = FF
        RETURN
        END

        real*8 FUNCTION F0_6125(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  1.086889851d0 ,
     2  1.076687383d0 ,
     3  1.077678499d0 ,
     4  1.057952088d0 ,
     5  1.015695632d0 ,
     6  0.9447325962d0 ,
     7  0.8255469517d0 ,
     8  0.6367476595d0 ,
     9  0.1778052515d0 ,
     1  -1.755568101d0 /

        FT1 = -23.23460003d0
        FTN = -9.464257141d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F0_6125 = FF
        RETURN
        END

        real*8 FUNCTION F0_535937(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  1.678704299d0 ,
     2  1.673504092d0 ,
     3  1.67621401d0 ,
     4  1.65987765d0 ,
     5  1.624792995d0 ,
     6  1.551303809d0 ,
     7  1.434398974d0 ,
     8  1.25958189d0 ,
     9  0.8310957793d0 ,
     1  -0.9929182642d0 /

        FT1 = -14.51615036d0
        FTN = -8.961012201d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        F0_535937 = FF
        RETURN
        END

        real*8 FUNCTION Fnn0_535937(TT)

        implicit none
        real*8 TT,FF,FT1,FTN
        real*8 F(1:10)
        real*8 T(1:10)
        real*8 splint2

        data T /
     1  0.001013380868d0 ,
     2  0.005095864126d0 ,
     3  0.01076356182d0 ,
     4  0.01863200735d0 ,
     5  0.02955574365d0 ,
     6  0.0447211298d0 ,
     7  0.06577518408d0 ,
     8  0.0950044564d0 ,
     9  0.1583407607d0 ,
     1  0.3800178256d0 /
        data F /
     1  1.175078385d0 ,
     2  1.177124107d0 ,
     3  1.184414294d0 ,
     4  1.184523224d0 ,
     5  1.198222105d0 ,
     6  1.177668254d0 ,
     7  1.149931288d0 ,
     8  1.168228305d0 ,
     9  1.220567678d0 ,
     1  1.393788479d0 /

        FT1 = 0.2411995169d0
        FTN = 0.8077635873d0
        FF = splint2(T,F,10,FT1,FTN,TT)
        Fnn0_535937 = FF
        RETURN
        END

        real*8 FUNCTION FALL(i,TT)
        implicit none
        integer i
        real*8 TT
        real*8 F3_58114
        real*8 F3_05139
        real*8 F2_6
        real*8 F2_3
        real*8 F2_2
        real*8 F2_1
        real*8 F2
        real*8 F1_86
        real*8 F1_75
        real*8 F1_6
        real*8 F1_5
        real*8 F1_4
        real*8 F1_3
        real*8 F1_2
        real*8 F1_1
        real*8 F0_95
        real*8 F0_8
        real*8 F0_7
        real*8 F0_6125
        real*8 F0_535937
        if (i.eq.1) FALL = F3_58114(TT)
        if (i.eq.2) FALL = F3_05139(TT)
        if (i.eq.3) FALL = F2_6(TT)
        if (i.eq.4) FALL = F2_3(TT)
        if (i.eq.5) FALL = F2_2(TT)
        if (i.eq.6) FALL = F2_1(TT)
        if (i.eq.7) FALL = F2(TT)
        if (i.eq.8) FALL = F1_86(TT)
        if (i.eq.9) FALL = F1_75(TT)
        if (i.eq.10) FALL = F1_6(TT)
        if (i.eq.11) FALL = F1_5(TT)
        if (i.eq.12) FALL = F1_4(TT)
        if (i.eq.13) FALL = F1_3(TT)
        if (i.eq.14) FALL = F1_2(TT)
        if (i.eq.15) FALL = F1_1(TT)
        if (i.eq.16) FALL = F0_95(TT)
        if (i.eq.17) FALL = F0_8(TT)
        if (i.eq.18) FALL = F0_7(TT)
        if (i.eq.19) FALL = F0_6125(TT)
        if (i.eq.20) FALL = F0_535937(TT)
        RETURN
        END

        real*8 FUNCTION FnnALL(i,TT)

        implicit none
        integer i
        real*8 TT
        real*8 Fnn3_58114
        real*8 Fnn0_535937
        if (i.eq.1) FnnALL = Fnn3_58114(TT)
        if (i.eq.20) FnnALL = Fnn0_535937(TT)
        RETURN
        END

        real*8 FUNCTION nnALL(i)
        implicit none
        integer i
        real*8 nn
        if (i.eq.1) nn= 0.005198117254d0
        if (i.eq.2) nn= 0.008402697116d0
        if (i.eq.3) nn= 0.01358286383d0
        if (i.eq.4) nn= 0.01962130473d0
        if (i.eq.5) nn= 0.02242039957d0
        if (i.eq.6) nn= 0.02577825447d0
        if (i.eq.7) nn= 0.02984155183d0
        if (i.eq.8) nn= 0.03709988454d0
        if (i.eq.9) nn= 0.04454482372d0
        if (i.eq.10) nn= 0.05828428092d0
        if (i.eq.11) nn= 0.07073553026d0
        if (i.eq.12) nn= 0.08700160883d0
        if (i.eq.13) nn= 0.1086629106d0
        if (i.eq.14) nn= 0.1381553325d0
        if (i.eq.15) nn= 0.1793631966d0
        if (i.eq.16) nn= 0.278445738d0
        if (i.eq.17) nn= 0.4662742473d0
        if (i.eq.18) nn= 0.6960128707d0
        if (i.eq.19) nn= 1.038946326d0
        if (i.eq.20) nn= 1.550846994d0
        nnALL = nn
        RETURN
        END

        integer FUNCTION nnn()
        implicit none
        nnn = 20
        RETURN
        END

C******************************************************************************
C*      Spline interface routines                                             *
C*                                                                            *
C*	Date: 06/26/06                                                        *
C*	Washington, DC                                                        *
C*	Jan Vorberger                                                         *   
C******************************************************************************

	REAL*8 FUNCTION splint2(x,y,n,yp1,ypn,xx) 
	INTEGER n,NMAX
	PARAMETER (NMAX=500)
	REAL*8 yp1,ypn,x(n),y(n),y2(NMAX),xx,yy
	call spline(x,y,n,yp1,ypn,y2)
	call splint(x,y,y2,n,xx,yy)
	splint2 = yy
	RETURN
	END

	REAL*8 FUNCTION splintnat(x,y,n,xx) 
	INTEGER n,NMAX
	PARAMETER (NMAX=500)
	REAL*8 yp1,ypn,x(n),y(n),y2(NMAX),xx,yy
	yp1 = 1.d30
	ypn = 1.d30
	call spline(x,y,n,yp1,ypn,y2)
	call splint(x,y,y2,n,xx,yy)
	splintnat = yy
	RETURN
	END

C*****Spline calculation*********************************************************
C*	x(1:n)...argument array
C*	y(1:n)...function array
C*	n...number of field points
C*	yp1,ypn...first derivatives at beginning and end of array
C*	y2(1:n)...2nd derivates of interpolating function (output)

	SUBROUTINE spline(x,y,n,yp1,ypn,y2)
	INTEGER n,NMAX
	REAL*8 yp1,ypn,x(n),y(n),y2(n)
	PARAMETER (NMAX=500)
	INTEGER i,k
	REAL*8 p,qn,sig,un,u(NMAX)
	if (yp1.gt..99e30) then
	   y2(1)=0.
	   u(1)=0.
	else
	   y2(1)=-0.5
	   u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
	endif
	do 11 i=2,n-1
	   sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
	   p=sig*y2(i-1)+2.
	   y2(i)=(sig-1.)/p
	   u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *  1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *  u(i-1))/p
 11	continue
	if (ypn.gt..99e30) then
	   qn=0.
	   un=0.
	else
	   qn=0.5
	   un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
	endif
	y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
	do 12 k=n-1,1,-1
	   y2(k)=y2(k)*y2(k+1)+u(k)
 12	continue
	return
	END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.

C*****Spline Interpolation*******************************************************
C*	xa(1:n)...argument array
C*	ya(1:n)...function array
C*	y2(1:n)...2nd derivates of interpolating function
C*	n...number of field points
C*	x...argument for desired function value
C*	y...(output) value at point x (by spline interpolation)

	SUBROUTINE splint(xa,ya,y2a,n,x,y)
	INTEGER n
	REAL*8 x,y,xa(n),y2a(n),ya(n),eps
	INTEGER k,khi,klo
	REAL*8 a,b,h

	eps = 0.1d0
	if (x.lt.xa(1)-eps*(xa(2)-xa(1))) then
	   write(*,*) " x-value below spline x-range ",xa(1),x,xa(n),n
	endif
	if (x.gt.xa(n)+eps*(xa(n)-xa(n-1))) then
	   write(*,*) " x-value above spline x-range ",xa(1),x,xa(n),n
	endif

	klo=1
	khi=n
 1	if (khi-klo.gt.1) then
	   k=(khi+klo)/2
	   if(xa(k).gt.x)then
	      khi=k
	   else
	      klo=k
	   endif
	   goto 1
	endif
	h=xa(khi)-xa(klo)
	if (h.eq.0.) then
	   write(*,*) 'bad xa input in splint'
	   stop
	endif
	a=(xa(khi)-x)/h
	b=(x-xa(klo))/h
	y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *  2)/6.
	return
	END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.
C*********************************************************************************
C*                                                                               *
C*       First principle equation of state of hydrogen-helium mixture            *
C*       Mixing ratio: 18 helium atoms per 220 hydrogen atoms                    *
C*                                                                               *
C*       Free energy interpolation code written by B. Militzer                   *
C*                                                                               *
C*       Copyright by B. Militzer                       Berkeley, CA 01-29-13    *
C*                                                                               *
C*********************************************************************************

      real*8 FUNCTION AUToK(T)
      implicit none
      real*8 T
      AUToK = 315774.660873d0 * T
      return
      end

      real*8 FUNCTION KToAU(TinK)
      implicit none
      real*8 TinK,AUToK
      KToAU = TinK / AUToK(1.0d0)
      return
      end
    
      real*8 FUNCTION RsFromDensity(nn)         
      implicit none
      real*8 nn,rs,pi
      pi = 3.14159265358979323846d0
      rs = (3.0d0/(4.0d0*pi*nn))**(1.0d0/3.0d0)
      RsFromDensity = rs
      return
      end

      real*8 FUNCTION DensityFromRs(rs)         
      implicit none
      real*8 nn,rs,pi
      pi = 3.14159265358979323846d0
      nn = 3.0d0/(4.0d0*pi*rs*rs*rs)
      DensityFromRs = nn
      return
      end

      real*8 FUNCTION DensityToGCC(nn)         
      implicit none
      real*8 nn
C* the mixing ratio (110+9) enters here
      DensityToGCC = 11.2935573673d0 * nn * (110d0+4d0*9d0)/128d0
      return
      end

      real*8 FUNCTION DensityFromGCC(nn)         
      implicit none
      real*8 DensityToGCC,nn
      DensityFromGCC = nn/DensityToGCC(1.0d0)
      return
      end

      real*8 FUNCTION PressureToGPa(P)         
      implicit none
      real*8 P
      PressureToGPa = 29421.009874d0*P
      return
      end

      real*8 FUNCTION EnergyToeV(E)         
      implicit none
      real *8 E
      EnergyToeV = E * 27.2113834279d0
      return
      end

C****************************************************************************
      real*8 FUNCTION FreeEnergy(nn,T)         
      implicit none
      
      integer i,imin,imax,n,NMAX,nnn
      PARAMETER (NMAX=500)
      real*8 nn,T,nns(NMAX),ffs(NMAX)
      real*8 rs,RsFromDensity,TinK
      real*8 AUToK
      real*8 eps,fp,fm
      real*8 nnALL,FALL,FnnALL,splint2
      real*8 fn1,fnn,f

      eps = 0.01d0
      fp  = 1.0d0+eps
      fm  = 1.0d0-eps

      rs = RsFromDensity(nn)
      TinK = AUToK(T)

      if (rs*fp .lt. RsFromDensity(nnALL(nnn()))) then
         write(*,*) "Density is too large: rs=",rs,
     x        RsFromDensity(nnALL(nnn()))
         stop
      endif

      if (rs*fm.gt.RsFromDensity(nnALL(1))) then
         write(*,*) "Density is too small: rs=",rs,
     x        RsFromDensity(nnALL(1))
         stop
      endif

      if (TinK*fp.lt.500d0) then
         write(*,*) "Tempature is too low",TinK
         stop
      endif

      if (TinK*fm.gt.120000d0) then
         write(*,*) "Tempature is too high",TinK
         stop
      endif

      imin = 1
      imax = nnn()

      n=0
      do i=imin,imax
         n=n+1
         nns(n) = nnALL(i)
         ffs(n) = FALL(i,T)
      enddo
      fn1 = FnnALL(imin,T)
      fnn = FnnALL(imax,T)
      f = splint2(nns,ffs,n,fn1,fnn,nn)

      if (nn*fp.lt.nns(1).or.nn*fm.gt.nns(n)) then
         write(*,*) "Spline problem:",TinK,RsFromDensity(nns(1)),
     *        RsFromDensity(nn),
     *        RsFromDensity(nns(n)),imin,imax
      endif

      FreeEnergy = f
      return
      end

C****************************************************************************

      real*8 FUNCTION Pressure(nn,T)
      implicit none
      real*8 nn,T,dnn,nn1,nn2,F1,F2,dFdnn,P,eps
      real*8 FreeEnergy
      eps = 1.0d-6
      dnn = nn*eps
      nn1 = nn-dnn
      nn2 = nn+dnn
      F1 = FreeEnergy(nn1,T)
      F2 = FreeEnergy(nn2,T)
      dFdnn = (F2-F1)/(2.0d0*dnn)
      P = dFdnn*nn*nn
      Pressure = P
      return
      end

      real*8 FUNCTION Entropy(nn,T)
      implicit none
      real*8 nn,T,dT,T1,T2,F1,F2,dFdT,FreeEnergy,eps
      eps = 1.0d-6
      dT = T*eps
      T1 = T-dT
      T2 = T+dT
      F1 = FreeEnergy(nn,T1)
      F2 = FreeEnergy(nn,T2)
      dFdT = (F2-F1)/(2.0d0*dT)
      Entropy = -dFdT
      return
      end

      real*8 FUNCTION InternalEnergy(nn,T)
      implicit none
      real*8 nn,T,E,FreeEnergy,Entropy
      E = FreeEnergy(nn,T) + T * Entropy(nn,T)
      InternalEnergy = E
      return
      end
      
C****************************************************************************
      real*8 FUNCTION MinimumTemperatureForDensity(nn)
      implicit none
      real*8 rs,RsFromDensity,nn
      real*8 T,TinK,AUToK,KToAU
      rs = RsFromDensity(nn)
      
      TinK = 500.0d0

      T = KToAU(TinK)
      MinimumTemperatureForDensity = T
      return
      end


      real*8 FUNCTION MaximumTemperatureForDensity(nn)
      implicit none
      real*8 rs,RsFromDensity,nn
      real*8 T,TinK,AUToK,KToAU
      rs = RsFromDensity(nn)
      
      TinK = 120000.0d0

      T = KToAU(TinK)
      MaximumTemperatureForDensity = T
      return
      end

      PROGRAM PrintEOS
      implicit none
      integer i,j,np
      real*8 nni,rs,rho
      real*8 T1,T2,T,TinK,AUToK,KToAU
      real*8 F,E,S,P
      real*8 PinGPa,EineV,FineV
      real*8 nnALL,RsFromDensity,DensityToGCC
      real*8 MinimumTemperatureForDensity,MaximumTemperatureForDensity
      real*8 FreeEnergy,Pressure,InternalEnergy,Entropy
      real*8 PressureToGPa,EnergyToeV

      write(*,*)  
      write(*,*) '#####################################################'
      write(*,*) '#                                                   #'
      write(*,*) '#   Equation of state of hydrogen-helium mixtures   #'
      write(*,*) '#                                                   #'
      write(*,*) '# Mixing ratio: 18 helium and 220 hydrogen atoms    #'
      write(*,*) '#               chosen to mimic Jupiter             #' 
      write(*,*) '# x = 2N_He/(2N_He+N_H) = 18/128 = 0.140625         #'
      write(*,*) '# Y =  M_He/( M_He+M_H) = 36/146 = 0.246575         #'
      write(*,*) '#                                                   #'
      write(*,*) '# Free energy interpolation code                    #'
      write(*,*) '# written by B. Militzer         Berkeley, 01-29-13 #'
      write(*,*) '#####################################################'

      do i=2,17 ! print from 2...(n-2) rather from 0...n to excluded two padded densities at low and high density
         nni = nnALL(i)
         rs  = RsFromDensity(nni)
         rho = DensityToGCC(nni)
         T1  = MinimumTemperatureForDensity(nni)
         T2  = MaximumTemperatureForDensity(nni)
         np = (T2-T1)/KToAU(500.0d0)+0.5d0
         write(*,*) 
         do j=0,np
c            T = T1*exp(log(T2/T1)*dble(j)/dble(np))
            T = T1+(T2-T1)*dble(j)/dble(np)
            TinK = AUToK(T)
            F = FreeEnergy(nni,T)
            P = Pressure(nni,T)
            E = InternalEnergy(nni,T)
            S = Entropy(nni,T)
            PinGPa = PressureToGPa(P)
            EineV  = EnergyToeV(E)
            FineV  = EnergyToeV(F)
c            write(*,*) rs,TinK,PinGPa,E,EineV
 10         format(a,f8.6,a,f8.6,a,f10.3,a,f10.6,a,f13.6,
     *           a,f12.6,a,f10.6)
            write(*,10) 
     *           'rs= ',rs,
     *           '  rhoHHe(g/cc)= ',rho,
     *           '  T(K)= ',TinK,
C    *           '  E(Ha/el)= ',E,
     *           '  E(eV/el)= ',EineV,
     *           '  P(GPa)= ',  PinGPa,
C    *           '  F(Ha/el)= ',F,
     *           '  F(eV/el)= ',FineV,
     *           '  S(kb/el)= ',S

         enddo
      enddo
      stop
      end
            
