00001
00002 include '../Module/phys_constant.f90'
00003 include '../Module/def_quantities.f90'
00004 include '../Module/def_quantities_derived.f90'
00005 include '../Subroutine/IO_input_physq_plot.f90'
00006 include '../Subroutine/IO_output_physq_plot.f90'
00007 include '../Analysis/Subroutine/IO_output_physq_paper.f90'
00008 include '../Analysis/Module/def_physq_to_array.f90'
00009
00010
00011
00012
00013
00014 PROGRAM find_secular_point
00015
00016 use def_physq_to_array
00017 implicit none
00018 integer :: flag_restmass = 0
00019 integer :: total_num_sol, nn = 100
00020 integer :: ii, jj, num_physq, iseq_end, iseq_sig
00021 real(8) :: katamuki, param
00022
00023 do ii = 1, nn
00024 jj = ii
00025 call IO_input_physq_plot(jj,flag_restmass)
00026 if (flag_restmass.eq.9999) then
00027 total_num_sol = ii - 1
00028 exit
00029 end if
00030 call physq_to_array(ii,flag_restmass,num_physq)
00031 end do
00032
00033 if (physq_array(9,1).gt.physq_array(9,2)) then
00034 write(6,*) ' Axis ratio decreasing order '
00035 iseq_end = 0
00036 iseq_sig = 1
00037 else
00038 write(6,*) ' Axis ratio increasing order '
00039 iseq_end = total_num_sol + 1
00040 iseq_sig = - 1
00041 end if
00042
00043 physq_array(:,iseq_end) = 0.0d0
00044 read(5,*) param
00045
00046 do jj = 3, num_physq
00047 katamuki = (physq_array(jj,iseq_end+iseq_sig*2) &
00048 & - physq_array(jj,iseq_end+iseq_sig*3)) &
00049 & /(physq_array( 9,iseq_end+iseq_sig*2) &
00050 & - physq_array( 9,iseq_end+iseq_sig*3))
00051 physq_array(jj,iseq_end) = katamuki &
00052
00053 & *(param - physq_array( 9,iseq_end+iseq_sig*1)) &
00054 & + physq_array(jj,iseq_end+iseq_sig*1)
00055 end do
00056
00057 call array_to_physq(iseq_end+iseq_sig*1,flag_restmass,num_physq)
00058 call IO_output_physq_plot(1,flag_restmass)
00059 call array_to_physq(iseq_end,flag_restmass,num_physq)
00060 call IO_output_physq_plot(2,flag_restmass)
00061
00062 call IO_output_physq_paper(1)
00063
00064 END PROGRAM find_secular_point
