1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
subroutine bisection(lr,lt0,y,eps0,lt,bivar,keyerrt,keybisect)
use eosmodule
implicit none
real*8 lr,lt0,y,eps0,lt
integer keyerrt
integer keybisect !this doesn't do anything...
! 1 -> operate on energy
! 2 -> operate on entropy
!temporary vars
real*8 lt1,lt2,ltmin,ltmax
real*8 f1,f2,fmid,dlt,ltmid
real*8 d1,d2,d3,tol
real*8 f1a,f2a
integer bcount,i,itmax,maxbcount
real*8 bivar(*)
bcount = 0
maxbcount = 80
tol=1.d-9 ! need to find energy to less than 1 in 10^-9
itmax=50
ltmax=logtemp(ntemp)
ltmin=logtemp(1)
lt = lt0
lt1 = dlog10(min(10.0d0**ltmax,1.10d0*(10.0d0**lt0)))
lt2 = dlog10(max(10.0d0**ltmin,0.90d0*(10.0d0**lt0)))
call findthis(lr,lt1,y,f1a,bivar,d1,d2,d3)
call findthis(lr,lt2,y,f2a,bivar,d1,d2,d3)
f1=f1a-eps0
f2=f2a-eps0
keyerrt=0
do while(f1*f2.ge.0.0d0)
bcount=bcount+1
lt1=dlog10(min(10.0d0**ltmax,1.2d0*(10.0d0**lt1)))
lt2=dlog10(max(10.0d0**ltmin,0.8d0*(10.0d0**lt2)))
call findthis(lr,lt1,y,f1a,bivar,d1,d2,d3)
call findthis(lr,lt2,y,f2a,bivar,d1,d2,d3)
f1=f1a-eps0
f2=f2a-eps0
if(bcount.ge.maxbcount) then
keyerrt=667
return
endif
enddo
if(f1.lt.0.0d0) then
lt=lt1
dlt=dlog10((10.0D0**lt2)-(10.0d0**lt1))
else
lt=lt2
dlt=dlog10((10.0D0**lt1)-(10.0d0**lt2))
endif
do i=1,itmax
dlt=dlog10((10.0d0**dlt)*0.5d0)
ltmid=dlog10(10.0d0**lt+10.0d0**dlt)
call findthis(lr,ltmid,y,f2a,bivar,d1,d2,d3)
fmid=f2a-eps0
if(fmid.le.0.0d0) lt=ltmid
if(abs(1.0d0-f2a/eps0).lt.tol) then
lt=ltmid
return
endif
enddo
end subroutine bisection
subroutine bisection_rho(lr0,lt,y,lpressin,lr,bivar,keyerrr,keybisect)
use eosmodule
implicit none
real*8 lr,lr0,y,lpressin,lt
integer keyerrr
integer keybisect !this doesn't do anything
! 3 -> operate on pressure
!temporary vars
real*8 lr1,lr2,lrmin,lrmax
real*8 f1,f2,fmid,dlr,lrmid
real*8 d1,d2,d3,tol
real*8 f1a,f2a
integer bcount,i,itmax,maxbcount
real*8 bivar(*)
bcount = 0
maxbcount = 80
tol=1.d-9 ! need to find energy to less than 1 in 10^-9
itmax=50
lrmax=logrho(nrho)
lrmin=logrho(1)
lr = lr0
lr1 = dlog10(min(10.0d0**lrmax,1.10d0*(10.0d0**lr0)))
lr2 = dlog10(max(10.0d0**lrmin,0.90d0*(10.0d0**lr0)))
call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3)
call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3)
f1=f1a-lpressin
f2=f2a-lpressin
keyerrr=0
do while(f1*f2.ge.0.0d0)
bcount=bcount+1
lr1=dlog10(min(10.0d0**lrmax,1.2d0*(10.0d0**lr1)))
lr2=dlog10(max(10.0d0**lrmin,0.8d0*(10.0d0**lr2)))
call findthis(lr1,lt,y,f1a,bivar,d1,d2,d3)
call findthis(lr2,lt,y,f2a,bivar,d1,d2,d3)
f1=f1a-lpressin
f2=f2a-lpressin
if(bcount.ge.maxbcount) then
keyerrr=667
return
endif
enddo
if(f1.lt.0.0d0) then
lr=lr1
dlr=dlog10((10.0D0**lr2)-(10.0d0**lr1))
else
lr=lr2
dlr=dlog10((10.0D0**lr1)-(10.0d0**lr2))
endif
do i=1,itmax
dlr=dlog10((10.0d0**dlr)*0.5d0)
lrmid=dlog10(10.0d0**lr+10.0d0**dlr)
call findthis(lrmid,lt,y,f2a,bivar,d1,d2,d3)
fmid=f2a-lpressin
if(fmid.le.0.0d0) lr=lrmid
if(abs(1.0d0-f2a/lpressin).lt.tol) then
lr=lrmid
return
endif
enddo
end subroutine bisection_rho
|
