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
|
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
! 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
|
