00001 subroutine calc_circular_orbit(total_iteration)
00002 use phys_constant, only : long
00003 use def_bh_parameter, only : ome_bh
00004 use def_quantities
00005 implicit none
00006 real(long) :: xa,xb,fa,fb,eps,delta,f1,f2
00007 integer :: total_iteration, iter_count, i, icycle, Nmax=50
00008 eps = 1.0d-07
00009 total_iteration = 0
00010
00011 xa = ome_bh
00012 call iteration_BBH_CF(iter_count)
00013 total_iteration = total_iteration + iter_count
00014 write(6,'(a21,i5)') '--- total iteration =', total_iteration
00015 call calc_physical_quantities_BBH_CF
00016 call save_solution(0)
00017 call write_last_physq
00018 call printout_physq_console_BBH
00019 fa = (admmass - komarmass)/komarmass
00020 ome_bh =ome_bh + 0.2d0*ome_bh
00021
00022 xb = ome_bh
00023 call iteration_BBH_CF(iter_count)
00024 total_iteration = total_iteration + iter_count
00025 write(6,'(a21,i5)') '--- total iteration =', total_iteration
00026 call calc_physical_quantities_BBH_CF
00027 call save_solution(1)
00028 call write_last_physq
00029 call printout_physq_console_BBH
00030 fb = (admmass - komarmass)/komarmass
00031
00032 write(6,*) '--------- Initial values for Omega and f(Omega)=(admmass-komarmass)/komarmass ---------------------------'
00033 write(6,'(a7,i2,1p,2e20.12)') 'cycle =',0,xa,fa
00034 write(6,'(a7,i2,1p,2e20.12)') 'cycle =',1,xb,fb
00035
00036 if ( dabs(fa)>dabs(fb) ) then
00037 f1 = xa
00038 xa = xb
00039 xb = f1
00040
00041 f1 = fa
00042 fa = fb
00043 fb = f1
00044 endif
00045
00046 icycle = 1
00047 do i = 2,Nmax
00048 icycle = icycle + 1
00049 if ( dabs(fa)>dabs(fb) ) then
00050 f1 = xa
00051 xa = xb
00052 xb = f1
00053 f1 = fa
00054 fa = fb
00055 fb = f1
00056 end if
00057 delta = (xb-xa)/(fb-fa)
00058 xb = xa
00059 fb = fa
00060 delta = delta*fa
00061 if ( dabs(delta) < eps ) then
00062 write(6,*) 'Convergence'
00063 write(6,*) '------------------------------------------------------------------'
00064 write(6,*) 'Omega orbit = ',ome_bh
00065 call printout_physq_console_BBH
00066 return
00067 end if
00068 if (i.lt.10) then
00069 xa = xa - i*delta/10.0d0
00070 else
00071 xa = xa - delta
00072 end if
00073 ome_bh = xa
00074 write(6,*) '------------------------------------------------------------------'
00075 write (6,'(a7,i2,a20,e20.12)') 'cycle =', icycle, ' New Omega = ',ome_bh
00076 call iteration_BBH_CF(iter_count)
00077 total_iteration = total_iteration + iter_count
00078 write(6,'(a21,i5)') '--- total iteration =', total_iteration
00079 call calc_physical_quantities_BBH_CF
00080 call save_solution(icycle)
00081 call write_last_physq
00082 fa = (admmass - komarmass)/komarmass
00083 write(6,'(a7,i2,1p,2e20.12)') 'cycle =',icycle,xa,fa
00084 call printout_physq_console_BBH
00085 end do
00086 write(6,*) 'No convergence after ',Nmax,' iterations.'
00087
00088 end subroutine calc_circular_orbit
