00001 subroutine calc_circular_orbit_BNS_mpt(total_iteration)
00002 use phys_constant, only : long, nmpt
00003 use def_quantities, only : admmass_asymp, komarmass_asymp
00004 implicit none
00005 real(long) :: xa,xb,fa,fb,eps,delta,f1,f2
00006 integer :: total_iteration, iter_count, i, icycle, Nmax=50
00007 eps = 1.0d-07
00008 total_iteration = 0
00009
00010 call iteration_BNS_CF_3mpt(iter_count)
00011 total_iteration = total_iteration + iter_count
00012 write(6,'(a21,i5)') '--- total iteration =', total_iteration
00013 call calc_physical_quantities_BNS_CF
00014 call save_solution_BNS_mpt(0)
00015 call write_last_physq_BNS_mpt
00016 call printout_physq_console_BNS_mpt
00017 call copy_def_matter_parameter_from_mpt(nmpt)
00018 xa = ome
00019 call copy_def_quantities_from_mpt(nmpt)
00020 fa = (admmass_asymp - komarmass_asymp)/komarmass_asymp
00021
00022 xb = xa + 0.2d0*xa
00023 call iteration_BNS_CF_3mpt(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_BNS_CF_3mpt(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_BNS_mpt
