Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/calc_eqsolver

00001 subroutine calc_eqsolver_TrpBH(rr,rini,rval)
00002   use phys_constant, only : long
00003   use def_bh_parameter, only : mass_pBH
00004   implicit none
00005   integer :: count, ii
00006   real(long) :: rr, rini, rval
00007   real(long) :: R1, R2, F1, F2, Rtmp, M, jacobian, fac0, error
00008   real(long) :: facp(5) = (/ 2.0d-1, 5.0d-1, 1.0d-0, 1.0d-0, 1.0d-0 /)
00009 !
00010 ! --- Solve parametric trunpet BH eqs using secant method.
00011 !
00012   M = mass_pBH
00013   R1 = rini
00014   R2 = rini*1.05d0
00015   F1 = fnc_R_trpBH(R1,M) - rr
00016 !
00017   count = 0
00018   do 
00019     count = count + 1
00020     ii = min0(5,count)
00021     fac0 = facp(ii)
00022 !
00023     F2 = fnc_R_trpBH(R2,M) - rr
00024     jacobian = (F2-F1)/(R2-R1)
00025     Rtmp = R2 - F2/jacobian
00026     R1 = R2
00027     F1 = F2
00028 !
00029     R2 = fac0*Rtmp + (1.0d0-fac0)*R2
00030     if (R2.le.3.0d0*M/2.0d0) R2 = 3.0d0*M/2.0d0
00031     error = 2.d0*dabs(R2 - R1)/(dabs(R2) + dabs(R1))
00032 ! 
00033     IF (error <= 1.d-14) EXIT 
00034     IF (count.ge.1000) exit
00035   end do
00036   IF (count.ge.1000) stop 'TrpBH failed'
00037 !
00038   rval   = R2
00039 !
00040   contains
00041     function fnc_R_trpBH(R,M)
00042       real(long) :: R, M, fnc_R_trpBH
00043       fnc_R_trpBH = 0.25d0*(2.0d0*R+M+dsqrt(4.0d0*R**2+4.0d0*M*R+3.0d0*M**2)) &
00044       &           *((4.0d0+3.0d0*dsqrt(2.0d0))*(2.0d0*R-3.0d0*M) &
00045       &  /(8.0d0*R+6.0d0*M+3.0d0*dsqrt(8.0d0*R**2+8.0d0*M*R+6.0d0*M**2)) &
00046       &            )**(1.0d0/dsqrt(2.0d0))
00047     end function fnc_R_trpBH
00048 end subroutine calc_eqsolver_TrpBH
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/calc_eqsolver_TrpBH.f90
