diff --git a/COSMOD_Catalogs.f95 b/COSMOD_Catalogs.f95 new file mode 100644 index 0000000..31f240c --- /dev/null +++ b/COSMOD_Catalogs.f95 @@ -0,0 +1,16 @@ +!=- fortran-libraries +!=- © Stanislav Shirokov, 2014-2020 + + module COSMOD_catalogs + use global + use COSMOD_paths + use math + use GNUplot + + character(len) :: Data_catalogs(100) = 'NaN' + + contains + + + + end module diff --git a/COSMOD_Citadel.f95 b/COSMOD_Citadel.f95 index 5e5d59d..ac10acb 100644 --- a/COSMOD_Citadel.f95 +++ b/COSMOD_Citadel.f95 @@ -8,60 +8,123 @@ module COSMOD_citadel use COSMOD_graphics use COSMOD_tables use COSMOD_scripts - - character(len) :: OS + use COSMOD_overleaf contains + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine commands + call preparation + do + call input_command(commandN) + select case (commandN(1)) + !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! + case('mn2021a') + call mn2021a - subroutine preparation + case('mn2020b') + call mn2020b - call define_system - if (operation_system==1) OS = 'Windows' - write(*,*) 'OS: ',trim(OS) - call paths - call gen_folders + case('mn2020a') + call mn2020a - write(*,*) 'COSMOD citadel is run..' + case('example') + call example - end subroutine + !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! + case('exit') + exit + case default + write(*,*) 'unknown command' - subroutine commands + end select - call preparation + call set_default + + enddo + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2021a + !=- If there is more one argument in command -=1 + if ( word_search_array( commandN , 'fig1' ) ) call mn2021a_fig1 + !if ( word_search_array( commandN , 'tab1' ) ) call mn2021a_plot_tab1 + !if ( word_search_array( commandN , 'o' ) ) call mn2021a_overleaf + !=- default command without additional arguments -=1 + if ( commandN(2)=='NaN' ) then + call mn2021a_fig1 + !call mn2021a_plot_tab1 + !call mn2021a_overleaf + end if - do ; write(*,*) 'Enter the command:' ; read (*,*) command + end subroutine - select case (command) + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2020b + !=- If there is more one argument in command -=1 + if ( word_search_array( commandN , 'fig1' ) ) call mn2020b_plot_fig1 + if ( word_search_array( commandN , 'tab1' ) ) call mn2020b_plot_tab1 + if ( word_search_array( commandN , 'o' ) ) call mn2020b_overleaf + !=- default command without additional arguments -=1 + if ( commandN(2)=='NaN' ) then + call mn2020b_plot_tab1 + call mn2020b_plot_fig1 + call mn2020b_overleaf + end if + + end subroutine - case ('overleaf') ; call make_overleaf_figures_folder - case ('o') ; call make_overleaf_figures_folder + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - case ('cosmod') - call COSMOD_model_set - case ('c') - call COSMOD_model_set + subroutine mn2020a + !=- If there is more one argument in command -=1 + if ( word_search_array( commandN , 'fig3' ) ) call mn2020a_plot_fig3 + if ( word_search_array( commandN , 'o' ) ) call mn2020a_overleaf + !=- default command without additional arguments -=1 + if ( commandN(2)=='NaN' ) then + call mn2020a_plot_fig3 + call mn2020a_overleaf + end if - case ('fig1') ; call MN_Letter_plot_fig1 + end subroutine - case ('tab1') ; call MN_Letter_plot_tab1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - case ('exit') ; exit + subroutine preparation + call define_system + call paths + write(*,*) 'COSMOD citadel has runned' + end subroutine - case default ; write(*,*) 'unknown command' - - end select + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - enddo + subroutine input_command(commandN) + character(len) command,commandN(:) + 33 continue ; write(*,*) '--' ; commandN(:)='NaN' + + read(*,'(A256)',err=33,end=44) command + read(command,*,iostat=iostat_value,err=33,end=44) commandN ; 44 continue end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine set_default !=- ? + + call MS_default + call overleaf_default + call math_default + + end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! end module +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/COSMOD_Config.f95 b/COSMOD_Config.f95 index c4ec0d1..9b3db19 100644 --- a/COSMOD_Config.f95 +++ b/COSMOD_Config.f95 @@ -1,8 +1,12 @@ !=- fortran-libraries !=- © Stanislav Shirokov, 2014-2020 -module COSMOD_config + module COSMOD_config + use global + + character(len),parameter :: CosMod_Version = '1.0' logical :: GRB_medians_flag = .true. -end module + end module +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/COSMOD_Graphics.f95 b/COSMOD_Graphics.f95 index a92b422..6835b07 100644 --- a/COSMOD_Graphics.f95 +++ b/COSMOD_Graphics.f95 @@ -12,14 +12,11 @@ module COSMOD_graphics character(len) caption - integer :: Omega_v_count , Omega_k_count , w_count - contains + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - - - subroutine MN_Letter_fig1 + subroutine mn2020b_fig1 integer i,j,k,ii,jj, col , start_position , i_max , k_max , j_max , n_LCDM , n_model, GNUfields_plots_i_start character(len) GRB_catalog_path , GRB_medians_path @@ -28,8 +25,8 @@ subroutine MN_Letter_fig1 call clear_plot_fields - pngterm = 'set term pngcairo font "Verdana,10" size 1920, 1080' !=- enhanced - epsterm = 'set term postscript enhanced color font "Verdana,16" size 8.5, 8.5' + pngterm = 'set term pngcairo font "Verdana,13" size 1920, 1080' !=- enhanced + epsterm = 'set term postscript enhanced color font "Verdana,13" size 8.5, 8.5' GNUfields(32) = 'set termoption dashed' @@ -41,23 +38,23 @@ subroutine MN_Letter_fig1 GNUplot_colors(2) = 'red' GNUplot_colors(3) = 'blue' - i_max = w_count - k_max = Omega_k_count - j_max = Omega_v_count + i_max = count_noNaN(MS_EoS_w) + k_max = count_noNaN(MS_Omega_k) + j_max = count_noNaN(MS_Omega_w) if ( MS_real_model_count /= i_max*j_max*k_max ) write(*,*) & 'COSMOD_Graphics: MN_Letter_fig1: incorrect model count' - n_LCDM = 8 + n_LCDM = 1 n_model = 0 do j=1,j_max !=- Omega_v 1-2*(j-1)*(j-3) do i=1,i_max !=- w do k=1,k_max !=- Omega_k & color - if (j==1) thickness = 2 - if (j==2) thickness = 4 - if (i==1) dash_type = 6 - if (i==2) dash_type = 1 + if (i==1) thickness = 4 + if (i==2) thickness = 2 + if (j==1) dash_type = 6 + if (j==2) dash_type = 1 !if (k==3) thickness = 2 !if (i==1 .and. j==2) dash_type = 6 !if (i==2 .and. j==2) dash_type = 1 @@ -117,17 +114,17 @@ subroutine MN_Letter_fig1 GNUfields(ylabel) = trim(MS_fields(1)) GNUfields(yrange) = '[38:*]' - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & + GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(MS_table_name)) // & '" u 1:' // trim( inttostr( col ) ) // ' w l ls ' // trim(inttostr(i)) case(2) - GNUfields(yrange) = '[*:*]' - GNUfields(yrange) = '[-1:1.6]' - if (GRB_medians_flag) GNUfields(yrange) = '[-2.3:1.6]' + !=- GNUfields(yrange) = '[*:*]' + GNUfields(yrange) = '[-1:1.4]' + if (GRB_medians_flag) GNUfields(yrange) = '[-2.3:1.4]' GNUfields(ylabel) = trim(MS_fields(2)) // ' = ' // & trim(MS_fields(3)) - GNUfields (plot1+(i-1)*2) = '"' // trim(slashfix(model_set_path)) // & + GNUfields (plot1+(i-1)*2) = '"' // trim(slashfix(MS_table_name)) // & '" u 1:(($' // trim( inttostr( col ) ) // ' - $' // & !'" u 1:(log10(10+$' // trim( inttostr( col ) ) // ' - $' // & !=- log10 trim(inttostr( 1+j+3*(n_LCDM-1) )) // ')) w l ls ' // trim(inttostr(i)) !=- log10 @@ -146,317 +143,91 @@ subroutine MN_Letter_fig1 caption='fig1a' endif graph_name = trim(MS_fields(4+1*(k-1))) // trim(caption) - GNUfields(extention_out_figure)='#eps' ; call plot(model_set_path) + GNUfields(extention_out_figure)='#eps' ; call plot(MS_table_name) graph_name = trim(MS_fields(4+1*(k-1))) // trim(caption) - GNUfields(extention_out_figure)='#png' ; call plot(model_set_path) - enddo - enddo - - end subroutine MN_Letter_fig1 - - - - - subroutine MN_Letter - integer i,j,k,ii,jj, col - call clear_plot_fields - - GNUfields(logscale) = '#set logscale' - GNUfields(format_y) = '#set format y "10^{%L}"' - GNUfields(xrange) = 'set xrange [0.1:*]' - - GNUfields(legend) = 'set key right bottom' - GNUfields(xlabel) = 'set xlabel "redshift z"' - GNUfields(grid) = 'set grid xtics ytics mxtics mytics' - GNUfields(mxtics) = 'set mxtics 4' - GNUfields(mytics) = 'set mytics 5' - - MS_fields(1) = '{/Symbol m}(z)' - !MS_fields(2) = 'log10 {/Symbol D}{/Symbol m}(z)' - !MS_fields(3) = 'log10 ( {/Symbol m}(z) - {/Symbol m}_{model-1}(z) )' - MS_fields(2) = '{/Symbol D}{/Symbol m}(z)' - MS_fields(3) = '{/Symbol m}(z) - {/Symbol m}_{model-1}(z)' - MS_fields(4) = 'mu(z)_' - MS_fields(5) = 'delta_mu(z)_' - MS_fields(6) = '{/Symbol D}{/Symbol m}(z)_' - do jj=1,3 - do ii=1,3 - MS_fields(7) = 'w_-1' - MS_fields(8) = 'w_-2' - MS_fields(9) = 'w_-3' - do k=1,2 - do j=3,3 - - GNUfields(title) = 'set title "'// trim(MS_fields(1+5*(k-1))) //' for different models"' - - do i=1,3 - col = 1+j+3*(i-1)+9*(ii-1)+27*(jj-1) - select case (k) - case (1) - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(1)) // '"' - GNUfields(yrange) = 'set yrange [38:58]' - - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:' // trim( inttostr( col ) ) // ' w l ls ' // trim(inttostr(i)) - GNUfields (title1+(i-1)*2) = ' title "model-' // trim(adjustl(inttostr(i))) // ':' // & - trim( model_set (i+3*(ii-1) + 9*(jj-1),2) ) // '"' - case(2) - GNUfields(yrange) = 'set yrange [0.97:1.08]' - !GNUfields(yrange) = 'set yrange [-0.05:-0.05]' - - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(2)) // ' = ' // & - trim(MS_fields(3)) // '"' - - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:(($' // trim( inttostr( col ) ) // '/$' // & !=- log10 - trim(inttostr( col - 3*(i-1) )) // ')) w l ls ' // trim(inttostr(i)) - GNUfields (title1+(i-1)*2) = ' title "model-' // trim(adjustl(inttostr(i))) // ':' // & - trim( model_set (i+3*(ii-1) + 9*(jj-1),2) ) // '"' - end select - end do - GNUfields (plot1) = GNUfields (plot1)(2:len) - - !graph_name = trim(MS_fields(4+1*(k-1))) // 'Ov_' // trim(MS_parameters(jj)) // '_' // & - ! trim(MS_fields(6+ii)) - ! GNUfields(extention_out_figure)='#eps' ; call plot(model_set_path) - graph_name = trim(MS_fields(4+1*(k-1))) // 'Ov_' // trim(MS_parameters(jj)) // '_' // & - trim(MS_fields(6+ii)) - GNUfields(extention_out_figure)='#png' ; call plot(model_set_path) - enddo - enddo - enddo - enddo - end subroutine MN_Letter - - - - subroutine MN_Letter_all - integer i,j,k,ii,jj, col - call clear_plot_fields - - pngterm = 'set term pngcairo font "Verdana,10" size 1920, 1080' !=- enhanced - epsterm = 'set term postscript enhanced color font "Verdana,14" size 19.2, 10.8' - - do i=1,3 - GNUfields( 33+1+9*(i-1) ) = 'set linestyle ' // trim(inttostr( 1+9*(i-1)) ) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "red"' - GNUfields(33+2+9*(i-1)) = 'set linestyle ' // trim(inttostr(2+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "blue"' - GNUfields(33+3+9*(i-1)) = 'set linestyle ' // trim(inttostr(3+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "dark-green"' - GNUfields(33+4+9*(i-1)) = 'set linestyle ' // trim(inttostr(4+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "yellow"' - GNUfields(33+5+9*(i-1)) = 'set linestyle ' // trim(inttostr(5+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "orange"' - GNUfields(33+6+9*(i-1)) = 'set linestyle ' // trim(inttostr(6+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "black"' - GNUfields(33+7+9*(i-1)) = 'set linestyle ' // trim(inttostr(7+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "purple"' - GNUfields(33+8+9*(i-1)) = 'set linestyle ' // trim(inttostr(8+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "brown"' - GNUfields(33+9+9*(i-1)) = 'set linestyle ' // trim(inttostr(9+9*(i-1))) // & - ' lw 1 pt 7 ps 0.7 dt ' // trim(inttostr(4+i)) // ' lc rgb "dark-cyan"' - end do - GNUfields(32) = 'set termoption dashed' - - GNUfields(logscale) = 'set logscale x' - GNUfields(format_y) = '#set format y "10^{%L}"' - GNUfields(xrange) = 'set xrange [0.1:*]' - - GNUfields(legend) = 'set key top left' - GNUfields(xlabel) = 'set xlabel "redshift z"' - GNUfields(grid) = 'set grid xtics ytics mxtics mytics' - GNUfields(mxtics) = 'set mxtics 10' - GNUfields(mytics) = 'set mytics 5' - - MS_fields(1) = '{/Symbol m}(z)' - !MS_fields(2) = 'log10 {/Symbol D}{/Symbol m}(z)' - !MS_fields(3) = 'log10 ( {/Symbol m}(z) - {/Symbol m}_{model-1}(z) )' - MS_fields(2) = '{/Symbol D}{/Symbol m}(z)' - MS_fields(3) = '{/Symbol m}(z) - {/Symbol m}_{model-1}(z)' - MS_fields(4) = 'mu(z)_' - MS_fields(5) = 'delta_mu(z)_' - MS_fields(6) = '{/Symbol D}{/Symbol m}(z)_' - - do k=1,2 - do j=3,3 - - GNUfields(title) = 'set title "'// trim(MS_fields(1+5*(k-1))) //' for different models"' - plot1 = 100 - title1 = 101 - - do i=1,27 - col = 1+j+3*(i-1) - select case (k) - case (1) - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(1)) // '"' - GNUfields(yrange) = 'set yrange [*:*]' - - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:' // trim( inttostr( col ) ) // ' w l ls ' // trim(inttostr(i)) - GNUfields (title1+(i-1)*2) = ' title "model-' // trim(adjustl(inttostr(i))) // ':' // & - trim( model_set ( i,2 )) // '"' - case(2) - GNUfields(yrange) = 'set yrange [*:*]' - !GNUfields(yrange) = 'set yrange [-0.05:-0.05]' - - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(2)) // ' = ' // & - trim(MS_fields(3)) // '"' - - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:(($' // trim( inttostr( col ) ) // ' - $' // & !=- log10 - trim(inttostr( 1+j )) // ')) w l ls ' // trim(inttostr(i)) - GNUfields (title1+(i-1)*2) = ' title "model-' // trim(adjustl(inttostr(i))) // ':' // & - trim( model_set (i,2) ) // '"' - end select - end do - GNUfields (plot1) = GNUfields (plot1)(2:len) - - graph_name = trim(MS_fields(4+1*(k-1))) // 'all_models' - GNUfields(extention_out_figure)='#eps' ; call plot(model_set_path) - graph_name = trim(MS_fields(4+1*(k-1))) // 'all_models' - GNUfields(extention_out_figure)='#png' ; call plot(model_set_path) + GNUfields(extention_out_figure)='#png' ; call plot(MS_table_name) enddo enddo - end subroutine MN_Letter_all + end subroutine mn2020b_fig1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - - subroutine MN_Letter_example - integer i,j,k,ii,jj, col , start_position , i_max , k_max , j_max , n_LCDM , n_model - character(len) GRB_catalog_path , GRB_medians_path - - GRB_catalog_path = 'C:\Users\Arhath\YandexDisk\Science\DATA\GRB_catalogs\GRBdata_Amati-193.dat' - GRB_medians_path = 'C:\Users\Arhath\YandexDisk\Science\DATA\GRB_catalogs\GRB_log_medians_data.dat' + subroutine plot_mn2021a_fig1(source_data) + character(len) source_data + character(2) :: delta_column = '7' call clear_plot_fields - pngterm = 'set term pngcairo font "Verdana,10" size 1920, 1080' !=- enhanced - epsterm = 'set term postscript enhanced color font "Verdana,16" size 8.5, 8.5' + graph_name = trim(adjustl(inttostr(approx_order))) // '_' // trim(name(source_data)) - GNUfields(32) = 'set termoption dashed' + pngterm = 'set term pngcairo enhanced font "Verdana,20" size 1400, 1050' + epsterm = 'set term postscript enhanced color font "Verdana,14" size 8.5, 6.3' - GNUplot_colors(1) = 'black' - GNUplot_colors(2) = 'red' - GNUplot_colors(3) = 'blue' - - start_position = 33 - i_max = w_count - k_max = Omega_k_count - j_max = Omega_v_count - - if ( MS_real_model_count /= i_max*j_max*k_max - 6 ) write(*,*) & - 'COSMOD_Graphics: MN_Letter_example: incorrect model count' + GNUfields (legend) = 'right bottom' - n_LCDM = 5 - - n_model = 0 - do j=1,j_max !=- Omega_v 1-2*(j-1)*(j-3) - if ( j /= 2 ) then - i_max = 1 - else - i_max = w_count - endif - do i=1,i_max !=- w - do k=1,k_max !=- Omega_k & color - if ( j==2 .and. i==1 ) then - thickness = 4 - else - thickness = 2 - endif - n_model = n_model + 1 - GNUfields( start_position+n_model ) = 'set linestyle ' // & - trim(inttostr( n_model )) // & - ' lw ' //trim(inttostr( thickness )) // ' pt 7 ps 0.7 dt ' // & - trim(inttostr( 2*(j-2)*(j-3)-(j-1)*(j-3)+6*(j-1)*(j-2)/2 )) // & - ' lc rgb "' // trim( GNUplot_colors(k) ) // '"' - enddo - enddo - enddo - GNUfields(start_position+n_model+1) = 'set linestyle ' // & - trim(inttostr( MS_real_model_count+1 )) // 'pt 7 ps 0.9 lc rgb "white"' - GNUfields(start_position+n_model+2) = 'set linestyle ' // & - trim(inttostr( MS_real_model_count+2 )) // 'pt 7 ps 0.9 lc rgb "dark-green"' + GNUfields (title) = ' approximation of the ' // trim(name(source_data)) // & + ' data (CosMod v1.0)' + GNUfields (xlabel) = 'z' + GNUfields (ylabel) = '{/Symbol m}' + GNUfields (logscale) = 'set logscale x' - GNUfields(logscale) = 'set logscale x' - GNUfields(format_y) = '#set format y "10^{%L}"' - GNUfields(xrange) = 'set xrange [0.1:10]' + GNUfields (ls1) = 'set linestyle 1 lw 1 pt 7 ps 0.7 lt rgb "blue"' + GNUfields (ls2) = 'set linestyle 2 lw 2 pt 7 ps 0.7 lt rgb "black"' + GNUfields (ls3) = 'set linestyle 3 lw 4 pt 7 ps 0.7 lt rgb "yellow"' + GNUfields (ls4) = 'set linestyle 4 lw 3 pt 7 ps 0.7 lt rgb "red"' - GNUfields(legend) = 'set key top left' - GNUfields(xlabel) = 'set xlabel "redshift z"' - GNUfields(grid) = 'set grid xtics ytics mxtics mytics' - GNUfields(mxtics) = 'set mxtics 10' - GNUfields(mytics) = 'set mytics 5' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:2:4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:5 w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:6 w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:7 w l ls 4' - MS_fields(1) = '{/Symbol m}(z)' + GNUfields (title1) = 'the ' // trim(name(source_data)) // ' data' + GNUfields (title2) = 'polylog approximation power ' // trim( adjustl( inttostr(approx_order) ) ) + GNUfields (title3) = 'first linear approximation (FLA)' + GNUfields (title4) = 'LCDM (70,0.7)' - MS_fields(2) = '{/Symbol D}{/Symbol m}(z)' - MS_fields(3) = '{/Symbol m}(z) - {/Symbol m}_{model-' // trim(inttostr(n_LCDM)) // '}(z)' - MS_fields(4) = 'mu(z)_' - MS_fields(5) = 'delta_mu(z)_' - MS_fields(6) = '{/Symbol D}{/Symbol m}(z)_' + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) - !MS_fields(2) = 'log_{10} ( 1 + {/Symbol D}{/Symbol m}(z) )' !=- log10 - !MS_fields(3) = 'log_{10} ( 1 + {/Symbol m}(z) - {/Symbol m}_{model-1}(z) )' !=- log10 - !MS_fields(6) = 'log_{10} 1+{/Symbol D}{/Symbol m}(z)_' !=- log10 + graph_name = 'delta_LCDM_' // trim(adjustl(inttostr(approx_order))) // '_' // trim(name(source_data)) - do k=1,2 !=- k=1 - absolutes, k=2 - residuals - do j=3,3 + GNUfields (ylabel) = '{/Symbol D}{/Symbol m} = {/Symbol m} - {/Symbol m}_{LCDM(70,0.7)}' - GNUfields(title) = 'set title "'// trim(MS_fields(1+5*(k-1))) //' for different models"' - plot1 = 100 - title1 = 101 + delta_column = '7' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:($2-$' // delta_column // '):4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:($5-$' // delta_column // ') w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:($6-$' // delta_column // ') w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:($7-$' // delta_column // ') w l ls 4' - do i=1,MS_real_model_count - col = 1+j+3*(i-1) + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) - select case (k) - case (1) - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(1)) // '"' - GNUfields(yrange) = 'set yrange [38:*]' + graph_name = 'delta_linear_' // trim(adjustl(inttostr(approx_order))) // '_' // trim(name(source_data)) - GNUfields (plot1+i*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:' // trim( inttostr( col ) ) // ' w l ls ' // trim(inttostr(i)) - GNUfields (plot1) = ', "' // trim(slashfix(GRB_catalog_path)) // & - '" u 2:3:4 w yerrorb ls ' // trim(inttostr( MS_real_model_count+1 )) - GNUfields (plot1+MS_real_model_count*2+2) = ', "' // trim(slashfix(GRB_medians_path)) // & - '" u 1:2:3 w yerrorb ls ' // trim(inttostr( MS_real_model_count+2 )) - case(2) - GNUfields(yrange) = 'set yrange [*:*]' - GNUfields(yrange) = 'set yrange [-2:4]' + GNUfields (ylabel) = '{/Symbol D}{/Symbol m} = {/Symbol m} - {/Symbol m}_{FLA}' - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(2)) // ' = ' // & - trim(MS_fields(3)) // '"' + delta_column = '6' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:($2-$' // delta_column // '):4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:($5-$' // delta_column // ') w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:($6-$' // delta_column // ') w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:($7-$' // delta_column // ') w l ls 4' - GNUfields (plot1+i*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:(($' // trim( inttostr( col ) ) // ' - $' // & - !'" u 1:(log10(10+$' // trim( inttostr( col ) ) // ' - $' // & !=- log10 - trim(inttostr( 1+j+3*(n_LCDM-1) )) // ')) w l ls ' // trim(inttostr(i)) - GNUfields (plot1) = ', "' // trim(slashfix(GRB_catalog_path)) // & - !'" u 2:5:4 w yerrorb ls ' // trim(inttostr( MS_real_model_count+1 )) - '" u 2:5 w p ls ' // trim(inttostr( MS_real_model_count+1 )) - !'" u 2:(log10(10+$5)) w p ls ' // trim(inttostr( MS_real_model_count+1 )) !=- log10 - GNUfields (plot1+MS_real_model_count*2+2) = ', "' // trim(slashfix(GRB_medians_path)) // & - '" u 1:5:3 w yerrorb ls ' // trim(inttostr( MS_real_model_count+2 )) - end select + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) - GNUfields (title1+i*2) = ' title "model-' // trim(adjustl(inttostr(i))) // ':' // & !=- - trim( model_set (i,2) ) // '"' - end do + graph_name = '5log_x_' // trim(adjustl(inttostr(approx_order))) // '_' // trim(name(source_data)) - GNUfields (plot1) = GNUfields (plot1)(2:len) + GNUfields (ylabel) = '{/Symbol D}{/Symbol m} = {/Symbol m} - {/Symbol m}_{55+5log x}' - GNUfields (title1) = ' notitle "GRBs(Amati-193)"' - GNUfields (title1+MS_real_model_count*2+2) = ' title "GRBs_{medians}(Amati-193)"' + delta_column = '8' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:($2-$' // delta_column // '):4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:($5-$' // delta_column // ') w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:($6-$' // delta_column // ') w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:($7-$' // delta_column // ') w l ls 4' - graph_name = trim(MS_fields(4+1*(k-1))) // 'all_models' - GNUfields(extention_out_figure)='#eps' ; call plot(model_set_path) - graph_name = trim(MS_fields(4+1*(k-1))) // 'all_models' - GNUfields(extention_out_figure)='#png' ; call plot(model_set_path) - enddo - enddo + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) + + end subroutine - end subroutine MN_Letter_example + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! end module COSMOD_graphics +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/COSMOD_Paths.f95 b/COSMOD_Paths.f95 index 380a192..3d60eb7 100644 --- a/COSMOD_Paths.f95 +++ b/COSMOD_Paths.f95 @@ -4,79 +4,34 @@ module COSMOD_paths use global - character(len) :: WorkDir='' , overleaf='' , MN_info , Models_folder , path_tabs + character(len) :: WorkDir , Models_folder , Tables_folder , overleaf_dir , & + + paper_name = 'NaN' contains + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine paths + !=- COSMOD folders WorkDir = trim( make_workdir() ) // 'Main_Workspace\' - Models_folder = trim(WorkDir) // 'Models\' - - overleaf = trim(WorkDir) // 'overleaf\MN_Letter\' - MN_info = trim(overleaf) // 'MN.info' - - path_tabs = trim(WorkDir) // 'Tables\' - - end subroutine - - - - subroutine gen_folders - - call system('mkdir ' // trim(Models_folder) // ' >> log.log ' ) - call system('mkdir ' // trim(path_tabs) // ' >> log.log ' ) - - call sleep (1) - - call system('cls') - - end subroutine gen_folders - - - - subroutine caption_MN_info - - open(unit_1,file=MN_info,status='replace') - write(unit_1,*) 'Anyone with this link can edit this project' - write(unit_1,*) 'https://www.overleaf.com/5957249232dfnmwjsjvdng' - write(unit_1,*) 'Anyone with this link can view this project' - write(unit_1,*) 'https://www.overleaf.com/read/ybwtxspzsdtb' - close(unit_1) - - end subroutine caption_MN_info - - - subroutine make_overleaf_figures_folder - - call system('echo off && MD ' // trim(overleaf) // ' >> log.log ' ) - - files( 1 ) = trim(WorkDir) // & - 'plots\MN_letter_model_set\MN_letter_model_set_mu(z)_all_models.eps' - new_files( 1 ) = 'Shirokov_fig1.eps' - - files( 2 ) = trim(WorkDir) // & - 'plots\MN_letter_model_set\MN_letter_model_set_delta_mu(z)_all_models.eps' - new_files( 2 ) = 'Shirokov_fig2.eps' - - files( 3 ) = trim(WorkDir) // & - 'plots\MN_letter_model_set\MN_letter_model_set_delta_mu(z)_fig1a.eps' - new_files( 3 ) = 'Shirokov_fig1a.eps' - - files( 4 ) = trim(WorkDir) // & - 'plots\MN_letter_model_set\MN_letter_model_set_delta_mu(z)_fig1b.eps' - new_files( 4 ) = 'Shirokov_fig1b.eps' + Models_folder = trim( WorkDir ) // 'Models\' + Tables_folder = trim( WorkDir ) // 'Tables\' - do i=1,N_files - if (files(i)/='') call system('copy '// trim(files(i)) //' '// trim(overleaf) // ' >> log.log ' ) - if (files(i)/='') call system('move '// trim(overleaf) // trim(name(files(i))) //'.eps '// & - trim(overleaf) // trim(new_files(i)) // ' >> log.log ' ) - enddo + overleaf_dir = trim( WorkDir ) // 'overleaf\' - call caption_MN_info + !=- creating folders from the Folders array + folders(1) = WorkDir + folders(2) = Models_folder + folders(3) = Tables_folder + folders(4) = overleaf_dir + call shell_MD_Tree end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! end module +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/COSMOD_Scripts.f95 b/COSMOD_Scripts.f95 index c5fdafd..47ac55e 100644 --- a/COSMOD_Scripts.f95 +++ b/COSMOD_Scripts.f95 @@ -1,191 +1,217 @@ !=- Fractal Dimension Estimation !=- © Stanislav Shirokov, 2014-2020 -module COSMOD_scripts - use math - use GNUplot - use cosmology + module COSMOD_scripts + use math + use GNUplot + use cosmology + + use COSMOD_config + use COSMOD_paths + use COSMOD_graphics + use COSMOD_tables + use COSMOD_catalogs + use COSMOD_overleaf - use COSMOD_config - use COSMOD_paths - use COSMOD_graphics + contains + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine example - contains + !=- models table + MS_table_name = trim(Models_folder) // 'example.dat' + !=- models parameters + MS_model(1) = 'wCDM' + MS_model(2) = 'FF' - subroutine MN_Letter_plot_fig1 + MS_Omega_k(1) = -0.1 + MS_Omega_k(2) = 0.0 + MS_Omega_k(3) = 0.1 - integer i,j,k,k_max,i_max,j_max , i_min , w , n_model - model_set_path = trim(WorkDir) // 'MN_letter_model_set.dat' + MS_EoS_w(1) = -1 + MS_EoS_w(2) = -2 - !=- write(*,*) ' the model_set format:' - !=- write(*,*) ' wCDM H_0 w Omega_v Omega_k' - !=- write(*,*) ' FF H_0' - !=- write(*,*) ' TL H_0' + MS_Omega_w(1) = 0.3 + MS_Omega_w(2) = 0.7 + + !=- calculating cosmological models + MS_recalculating = .true. + call make_model_set + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2021a_fig1 + integer i + + Data_catalogs(1) = 'C:\Users\Arhath\YandexDisk\Science\DATA\SN_catalogs\Pantheon.dat' + + columns_for_approx = '6 3 5' + approx_order = 8 + call approx_polylog(Data_catalogs(1)) + call plot_mn2021a_fig1(Data_catalogs(1)) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2020b_overleaf + + paper_name = 'mn2020b' + + overleaf_CanEdit = 'https://www.overleaf.com/5957249232dfnmwjsjvdng' + overleaf_ReadOnly = 'https://www.overleaf.com/read/ybwtxspzsdtb' + + files( 1 ) = 'Models\plots\mn2020b_fig1\mn2020b_fig1_delta_mu(z)_fig1a.eps' + new_files( 1 ) = 'Shirokov_fig1a.eps' + + files( 2 ) = 'Models\plots\mn2020b_fig1\mn2020b_fig1_delta_mu(z)_fig1b.eps' + new_files( 2 ) = 'Shirokov_fig1b.eps' + + call overleaf_making + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2020b_plot_tab1 + + MS_model(1) = 'wCDM' + MS_Hubble_0(1) = 70 + + MS_Omega_w(1) = 0.3d0 + MS_Omega_w(2) = 0.5d0 + MS_Omega_w(3) = 0.7d0 + MS_Omega_w(4) = 0.9d0 + + MS_EoS_w(1) = -2 + + MS_Omega_k(1) = -0.1d0 + MS_Omega_k(2) = 0.0d0 + MS_Omega_k(3) = 0.1d0 + + tab1_z(1 ) = 0.1d0 + tab1_z(2 ) = 0.2d0 + tab1_z(3 ) = 0.3d0 + tab1_z(4 ) = 0.5d0 + tab1_z(5 ) = 0.6d0 + tab1_z(6 ) = 0.8d0 + tab1_z(7 ) = 1.d0 + tab1_z(8 ) = 2.d0 + tab1_z(9 ) = 3.d0 + tab1_z(10) = 4.d0 + tab1_z(11) = 5.d0 + tab1_z(12) = 6.d0 + tab1_z(13) = 8.d0 + tab1_z(14) = 10.d0 + + call mn2020b_tab1 - MS_parameters(1)='0.3' - MS_parameters(2)='0.7' - - w_count = 2 - Omega_v_count = 2 - Omega_k_count = 3 - MS_real_model_count = w_count*Omega_k_count*Omega_v_count - - - n_model=0 - do j=1,Omega_v_count - do i=1,w_count - w=-i - do k=1,Omega_k_count - MS_O_k = -0.2 + 0.1*k - mm=4 - if (k>1) mm=5 - read(MS_parameters(j),*) MS_O_v - MS_O_m_bug = abs ( 1 + MS_O_k - MS_O_v ) !=- the Friedmann equation - - n_model = n_model + 1 - model_set ( n_model , 1 ) = & - ' wCDM 70 ' // trim(inttostr(w)) // ' ' // & - trim(MS_parameters(j)) // ' ' // trim(realtostrf(MS_O_k,4,1)) - model_set ( n_model , 2 ) = & - ' wCDM (w = ' // trim(inttostr(w)) // & - ', {/Symbol W}_{DE}=' // trim(MS_parameters(j)) // & - ', {/Symbol W}_m=' // trim(realtostrf(MS_O_m_bug,4,1)) // & - ', {/Symbol W}_k=' // trim(realtostrf(MS_O_k,mm,1)) // ') ' - enddo - enddo - enddo - - !=- the hand input: - !model_set - - !=- write(*,'(A)') model_set(:,1) ; pause - - call compute_model_set - call write_model_set + end subroutine mn2020b_plot_tab1 + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2020b_plot_fig1 + + !=- models table + MS_table_name = trim(Models_folder) // 'mn2020b_fig1.dat' + + !=- the second way setting Model_Set parameters + !=- models parameters + MS_model(1) = 'wCDM' + + MS_Omega_k(1) = 0.0 !=- first is LCDM + MS_Omega_k(2) = -0.1 + MS_Omega_k(3) = 0.1 + + MS_EoS_w(1) = -1 !=- first is LCDM + MS_EoS_w(2) = -2 + + MS_Omega_w(1) = 0.7 !=- first is LCDM + MS_Omega_w(2) = 0.3 + + !=- calculating cosmological models + !=- MS_recalculating = .true. + call make_model_set + + !=- plotting figures GRB_medians_flag = .false. - call MN_Letter_fig1 + call mn2020b_fig1 GRB_medians_flag = .true. - call MN_Letter_fig1 - !=- call MN_Letter_example + call mn2020b_fig1 + end subroutine - end subroutine MN_Letter_plot_fig1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine mn2020a_overleaf - subroutine COSMOD_model_set - integer i,j,k,k_max,i_max,j_max , i_min , w , n_model - model_set_path = trim(WorkDir) // 'MN_letter_model_set.dat' + paper_name = 'mn2020a' - !=- write(*,*) ' the model_set format:' - !=- write(*,*) ' wCDM H_0 w Omega_v Omega_k' - !=- write(*,*) ' FF H_0' - !=- write(*,*) ' TL H_0' + overleaf_CanEdit = 'https://www.overleaf.com/9181497425tkqhjdzzmbgg' + overleaf_ReadOnly = 'https://www.overleaf.com/read/jjrmkfbhhftv' + + files( 1 ) = 'Models\plots\mn2020a_fig3\mn2020a_fig3_r(z)-logscale.eps' + new_files( 1 ) = 'Shirokov_fig3a.eps' + + files( 2 ) = 'Models\plots\mn2020a_fig3\mn2020a_fig3_delta_r(z).eps' + new_files( 2 ) = 'Shirokov_fig3b.eps' - MS_parameters(1)='0.5' - MS_parameters(2)='0.7' - MS_parameters(3)='0.9' - - w_count = 2 - Omega_v_count = 3 - Omega_k_count = 3 - MS_real_model_count = w_count*Omega_k_count*Omega_v_count - 6 - - - n_model=0 - do j=1,Omega_v_count - if ( j /= 2 ) then - i_max = 1 - else - i_max = w_count - endif - do i=1,i_max - w=-i - if ( j /= 2 ) w=-2 - do k=1,Omega_k_count - MS_O_k = -0.2 + 0.1*k - mm=4 - if (k>1) mm=5 - read(MS_parameters(j),*) MS_O_v - MS_O_m_bug = abs ( 1 + MS_O_k - MS_O_v ) !=- the Friedmann equation - !write(*,*) MS_O_m_bug ; pause - n_model = n_model + 1 - model_set ( n_model , 1 ) = & - ' wCDM 70 ' // trim(inttostr(w)) // ' ' // & - trim(MS_parameters(j)) // ' ' // trim(realtostrf(MS_O_k,4,1)) - model_set ( n_model , 2 ) = & - ' wCDM (w = ' // trim(inttostr(w)) // & - ', {/Symbol W}_{DE}=' // trim(MS_parameters(j)) // & - ', {/Symbol W}_m=' // trim(realtostrf(MS_O_m_bug,4,1)) // & - ', {/Symbol W}_k=' // trim(realtostrf(MS_O_k,mm,1)) // ') ' - enddo - enddo - enddo - - !=- the hand input: - !model_set - - !=- write(*,'(A)') model_set(:,1) ; pause - - call compute_model_set - call write_model_set - call MN_Letter_example - - end subroutine COSMOD_model_set - - - - subroutine COSMOD_model_set_all - integer i,j,k - model_set_path = trim(WorkDir) // 'MN_letter_' // model_set_path + files( 3 ) = 'Models\plots\mn2020a_fig3\mn2020a_fig3_d_L(z)-logscale.eps' + new_files( 3 ) = 'Shirokov_fig3c.eps' + + files( 4 ) = 'Models\plots\mn2020a_fig3\mn2020a_fig3_delta_d_L(z).eps' + new_files( 4 ) = 'Shirokov_fig3d.eps' + + files( 5 ) = 'Models\plots\mn2020a_fig3\mn2020a_fig3_mu(z)-logscale.eps' + new_files( 5 ) = 'Shirokov_fig4a.eps' + + files( 6 ) = 'Models\plots\mn2020a_fig3\mn2020a_fig3_delta_mu(z).eps' + new_files( 6 ) = 'Shirokov_fig4b.eps' + + call overleaf_making + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2020a_plot_fig3 + + !=- models table + MS_table_name = trim(Models_folder) // 'mn2020a_fig3.dat' !=- write(*,*) ' the model_set format:' !=- write(*,*) ' wCDM H_0 w Omega_v Omega_k' !=- write(*,*) ' FF H_0' !=- write(*,*) ' TL H_0' - MS_parameters(1)='0.7' - MS_parameters(2)='0.5' - MS_parameters(3)='0.9' - - do k=1,3 - read(MS_parameters(k),*) MS_O_v - do i=1,3 - MS_O_m_bug = abs ( 1 + 0.0 - MS_O_v ) !=- the Friedmann equation - !write(*,*) MS_O_m_bug ; pause - model_set ( 1+3*(i-1)+9*(k-1) , 1 ) = ' wCDM 70 ' // trim(inttostr(-i)) // ' ' // & - trim(MS_parameters(k)) // ' 0.0 ' - model_set ( 1+3*(i-1)+9*(k-1) , 2 ) = ' the wCDM (w = ' // trim(inttostr(-i)) // & - ', {/Symbol W}_{/Symbol L}=' // trim(MS_parameters(k)) // & - ', {/Symbol W}_m=' // trim(realtostrf(MS_O_m_bug,3,1)) // & - ', {/Symbol W}_k= 0.0) model ' - - MS_O_m_bug = abs ( 1 + 0.1 - MS_O_v ) !=- the Friedmann equation - model_set ( 2+3*(i-1)+9*(k-1) , 1 ) = ' wCDM 70 ' // trim(inttostr(-i)) // ' ' // & - trim(MS_parameters(k)) // ' 0.1 ' - model_set ( 2+3*(i-1)+9*(k-1) , 2 ) = ' the wCDM (w = ' // trim(inttostr(-i)) // & - ', {/Symbol W}_{/Symbol L}=' // trim(MS_parameters(k)) // & - ', {/Symbol W}_m=' // trim(realtostrf(MS_O_m_bug,3,1)) // & - ', {/Symbol W}_k= 0.1) model ' - - MS_O_m_bug = abs ( 1 - 0.1 - MS_O_v ) !=- the Friedmann equation - model_set ( 3+3*(i-1)+9*(k-1) , 1 ) = ' wCDM 70 ' // trim(inttostr(-i)) // ' ' // & - trim(MS_parameters(k)) // ' -0.1 ' - model_set ( 3+3*(i-1)+9*(k-1) , 2 ) = ' the wCDM (w = ' // trim(inttostr(-i)) // & - ', {/Symbol W}_{/Symbol L}=' // trim(MS_parameters(k)) // & - ', {/Symbol W}_m=' // trim(realtostrf(MS_O_m_bug,3,1)) // & - ', {/Symbol W}_k=-0.1) model ' - enddo - enddo - - call compute_model_set - call write_model_set - call MN_Letter_all - - end subroutine COSMOD_model_set_all + !=- the first way setting Model_Set parameters + model_set ( 1 , 1 ) = ' wCDM 70 -1 0.7 0.0 ' + model_set ( 1 , 2 ) = ' {/Symbol L}CDM( {/Symbol W}_{/Symbol L} = 0.7 ) ' + + model_set ( 2 , 1 ) = ' wCDM 70 -1 1.0 0.0 ' + model_set ( 2 , 2 ) = ' {/Symbol L}CDM( {/Symbol W}_{/Symbol L} = 1.0 ), CSS ' + + model_set ( 3 , 1 ) = ' TL 70 ' + model_set ( 3 , 2 ) = ' the TL model ' + + model_set ( 4 , 1 ) = ' FF 70 ' + model_set ( 4 , 2 ) = ' the FF model ' + + model_set ( 5 , 1 ) = ' wCDM 70 -1 0.9 0.0 ' + model_set ( 5 , 2 ) = ' {/Symbol L}CDM( {/Symbol W}_{/Symbol L} = 0.9 ) ' + + model_set ( 6 , 1 ) = ' wCDM 70 -0.5 0.7 0.2 ' + model_set ( 6 , 2 ) = ' wCDM( w = -0.5, {/Symbol W}_{/Symbol L} = 0.7, {/Symbol W}_k = 0.2 ) ' + + call make_model_set + end subroutine mn2020a_plot_fig3 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! end module COSMOD_scripts +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/COSMOD_Tables.f95 b/COSMOD_Tables.f95 index 9ca0e21..2595a07 100644 --- a/COSMOD_Tables.f95 +++ b/COSMOD_Tables.f95 @@ -1,136 +1,85 @@ -!=- fortran-libraries -!=- © Stanislav Shirokov, 2014-2020 - - module COSMOD_tables - use global - use GNUplot - use cosmology - use math - - use COSMOD_config - use COSMOD_paths - - character(len) :: path_tab1_dat , path_tab1_tex - - integer,parameter :: range_z = 14 , range_model = 15 , range_omega_v = 5 , range_w = 1 , range_omega_k = 3 - - real(8) :: tab1_mu(range_z,1+range_model) , & - value_omega_v(range_omega_v) , & - value_w(range_w) , & - value_omega_k(range_omega_k) , & - - mu_LCDM - - contains - - - - subroutine MN_Letter_plot_tab1 - integer w,v,k - - path_tab1_dat = trim(path_tabs) // 'MN_Letter_tab1.dat' - path_tab1_tex = trim(path_tabs) // 'MN_Letter_tab1.tex' - - value_omega_v(1) = 0.9d0 - value_omega_v(2) = 0.7d0 - value_omega_v(3) = 0.5d0 - value_omega_v(4) = 0.3d0 - value_omega_v(5) = 0.1d0 - - value_w(1) = -2 - - value_omega_k(1) = -0.1d0 - value_omega_k(2) = 0.0d0 - value_omega_k(3) = 0.1d0 - - tab1_mu(1 ,1) = 0.1d0 - tab1_mu(2 ,1) = 0.2d0 - tab1_mu(3 ,1) = 0.3d0 - tab1_mu(4 ,1) = 0.5d0 - tab1_mu(5 ,1) = 0.6d0 - tab1_mu(6 ,1) = 0.8d0 - tab1_mu(7 ,1) = 1.d0 - tab1_mu(8 ,1) = 2.d0 - tab1_mu(9 ,1) = 3.d0 - tab1_mu(10,1) = 4.d0 - tab1_mu(11,1) = 5.d0 - tab1_mu(12,1) = 6.d0 - tab1_mu(13,1) = 8.d0 - tab1_mu(14,1) = 10.d0 - - do i=1,range_z - mu_LCDM = dm_wCDM( tab1_mu(i,1) , 7d1 , -1d0 , 0.3d0 , 0d0 ) - do w=1,range_w - do v=1,range_omega_v - do k=1,range_omega_k - - j = 1 + k + range_omega_k*(v-1) + range_omega_k*range_omega_v*(w-1) - - tab1_mu(i,j) = dm_wCDM( tab1_mu(i,1) , 7d1 , value_w(w) , & - 1 + value_omega_k(k) - value_omega_v(v) , value_omega_k(k)) - mu_LCDM - write(*,*) tab1_mu(i,1) , mu_LCDM , tab1_mu(i,j) , value_w(w) , & - 1 + value_omega_k(k) - value_omega_v(v) , value_omega_k(k) , & - dm_wCDM( tab1_mu(i,1) , 7d1 , value_w(w) , & - 1 + value_omega_k(k) - value_omega_v(v) , value_omega_k(k)) - end do - end do - end do - end do - - theformat = '(' // trim(inttostr(1+range_model)) // '(F5.2,1x))' - call write_tab1(path_tab1_dat) - - theformat = '(' // trim(inttostr(range_model)) // '(F5.2,1x,"&",1x),F5.2,1x,"\\")' - call write_tab1(path_tab1_tex) - - !call write_tab1b - - end subroutine - - - - subroutine write_tab1(path_table) - character(len) path_table ; integer i - - unit_1 = random_unit() - - open(unit_1,file=path_table,status='replace') - if (theformat=='') theformat='*' - write(unit_1,theformat) ( tab1_mu(i,:) , i=1 , range_z) - - close(unit_1) - - end subroutine - - - - subroutine write_tab1b - integer w,v,k - - theformat = '(i2,1x,"&",1x,F3.1,1x,"&",1x,F4.1,' // trim(inttostr(range_z)) // & - '(1x,"&",1x,F5.2),1x,"\\")' - - unit_1 = random_unit() - - open(unit_1,file=path_tab1_tex,status='replace') - if (theformat=='') theformat='*' - - do w=1,range_w - do v=1,range_omega_v - do k=1,range_omega_k - - j = 1 + k + range_omega_k*(v-1) + range_omega_k*range_omega_v*(w-1) - write(unit_1,theformat) int(value_w(w)) , value_omega_v(v) , value_omega_k(k) , & - tab1_mu(:,j) +!=- fortran-libraries +!=- © Stanislav Shirokov, 2014-2020 + + module COSMOD_tables + use global + use GNUplot + use cosmology + use math + + use COSMOD_config + use COSMOD_paths + + integer,parameter :: range_z = 1d2 + + character(len) :: path_tab1_dat , path_tab1_tex + + real(8) :: tab1_z(range_z) = NaN + + real(8),allocatable,dimension(:,:) :: tab1_mu + + contains + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine mn2020b_tab1 + integer z,w,v,k, model , count_z + real(8) mu_LCDM + + call cheking_MS_parameters + MS_real_model_count = MS_count_w * MS_count_O_w * MS_count_O_k + + count_z = count_noNaN(tab1_z) + + allocate(tab1_mu(count_z,1+MS_real_model_count)) + + path_tab1_dat = trim(Tables_folder) // 'mn2020b_tab1.dat' + path_tab1_tex = trim(Tables_folder) // 'mn2020b_tab1.tex' + + tab1_mu(:,1) = tab1_z(:) + do z=1,count_z + + mu_LCDM = dm_wCDM( tab1_mu(z,1) , 7d1 , -1d0 , 0.3d0 , 0d0 ) + + do w=1,MS_count_w + do v=1,MS_count_O_w + do k=1,MS_count_O_k + + model = 1 + k + MS_count_O_k*(v-1) + MS_count_O_k*MS_count_O_w*(w-1) + + tab1_mu( z , model ) = dm_wCDM( tab1_mu(z,1) , 7d1 , MS_EoS_w(w) , & + 1 + MS_Omega_k(k) - MS_Omega_w(v) , MS_Omega_k(k)) - mu_LCDM end do end do - end do - - close(unit_1) - - end subroutine - - - - end module + end do + end do + + theformat = '(' // trim(inttostr(1+MS_real_model_count)) // '(F5.2,1x))' + call write_tab1(path_tab1_dat) + + theformat = '(' // trim(inttostr(MS_real_model_count)) // '(F5.2,1x,"&",1x),F5.2,1x,"\\")' + call write_tab1(path_tab1_tex) + + deallocate(tab1_mu) + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine write_tab1(path_table) + character(len) path_table ; integer i + + unit_1 = random_unit() + + open(unit_1,file=path_table,status='replace') + if (theformat=='') theformat='*' + write(unit_1,theformat) ( tab1_mu(i,:) , i=1 , size(tab1_mu(:,1)) ) + + close(unit_1) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + end module +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/COSMOD_overleaf.f95 b/COSMOD_overleaf.f95 new file mode 100644 index 0000000..acb4abb --- /dev/null +++ b/COSMOD_overleaf.f95 @@ -0,0 +1,60 @@ +!=- fortran-libraries +!=- © Stanislav Shirokov, 2014-2020 + +module COSMOD_overleaf + use global + use COSMOD_paths + + character(len) :: overleaf_info , overleaf_CanEdit , overleaf_ReadOnly + + contains + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine overleaf_making + + character(len) filename, new_filename + + overleaf_info = trim(overleaf_dir) // trim(paper_name) // '\' + call create_folder(overleaf_info) + + forall ( i=1:N_files , files(i)/='NaN' ) files(i) = trim(WorkDir) // trim(files(i)) + + do i=1,N_files + + if (files(i)/='NaN') call copy_file( files(i) , overleaf_info ) + + filename = trim(overleaf_info) // trim(name(files(i))) //'.eps' + new_filename = trim(overleaf_info) // trim( new_files(i)) + if (files(i)/='NaN') call move_file( filename , new_filename ) + + enddo + + overleaf_info = trim(overleaf_info) // trim(paper_name) // '.info' + + unit_1 = random_unit() + open(unit_1,file=overleaf_info,status='replace') + write(unit_1,*) 'Anyone with this link can edit this project' + write(unit_1,*) trim(overleaf_CanEdit) + write(unit_1,*) 'Anyone with this link can view this project' + write(unit_1,*) trim(overleaf_ReadOnly) + write(unit_1,*) 'Local filepath: ', trim(overleaf_info) + close(unit_1) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine overleaf_default + + paper_name = 'NaN' + new_files(:) = 'NaN' + files (:) = 'NaN' + folders (:) = 'NaN' + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + end module +!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/Cosmology.f95 b/Cosmology.f95 index 2c4e315..52239ab 100644 --- a/Cosmology.f95 +++ b/Cosmology.f95 @@ -1,102 +1,231 @@ !=- fortran-libraries -!=- © Stanislav Shirokov, 2014-2020 +!=- © Stanislav Shirokov, 2014-2020 -!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 module cosmology use global use GNUplot - logical :: RDR_calculating = .true. , & !=- .false. - RDR_error = .true. + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + !=- cosmological constants + real(8),parameter :: c = 2.99792458d10 , & !=- cm/s + + L_sun = 3.827d33 , & !=- erg/s + M_sun = 4.83d0 , & + + O_3 = 0.3d0 , & + O_7 = 0.7d0 , & + H70 = 7d1 , & !=- km/s/Mpc + dl = 4421.71767d0 !=- c/H_0/1d5 in GCS + + !=- cosmological variables + real(8) :: z_max = 2d1 , & !=- redshift calculated maximally - integer,parameter :: N_models = 9 , & - N_RDR_grid = 2d3 , & !=- graph grid - N_RDR_grid_mode = 3 , & !=- progression power + H_0 = 70d0 , & !=- km/s/Mpc + O_m = 0.308d0 , & + O_w = 1d0 , & + O_k = 0d0 , & + q_0 = 0.5d0 - MS_grid = 2d2 , & - MS_model_count = 1d2 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - integer :: N_grid = 2d2 , & + !=- others variables + real(8) :: dz, x_max, I1 - RDR_type = 1 , & !=- 0 is LCDM(0.3,70) , 1 is LCDM(\O_m,H_0) , 2 is wCDM(w,\O_m,H_0,O_k) - MS_real_model_count + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - real(8) :: c = 2.99792458d10 , & !=- cm/s - H_0 = 70d0 , & !=- km/s/Mpc - H70 = 7d1 , & - dl = 4421.71767d0 , & !=- c/H_0/1d5 in GCS + !=- the redshift_distance_relation procedure + logical :: RDR_calculating = .true. , & !=- .false. + RDR_error = .true. - L_sun = 3.827d33 , & !=- erg/s - M_sun = 4.83d0 , & + integer,parameter :: N_RDR_grid = 2d3 , & !=- graph grid + N_RDR_grid_mode = 3 !=- progression power - q_0 = 0.5d0 , & - O_m = 0.308d0 , & - O_w = 1d0 , & - O_k = 0d0 , & - z_max = 2d1 , & - MS_z_max = 2d1 , & - RDR_z_max = 2d1 , & + real(8) :: RDR(2,N_RDR_grid) , & + RDR_z_max = 2d1 - O_3 = 0.3d0 , & - O_7 = 0.7d0 , & + integer :: RDR_type = 1 !=- 0 is LCDM(0.3,70) , 1 is LCDM(\O_m,H_0) , 2 is wCDM(w,\O_m,H_0,O_k) - MS_H_0 , MS_O_v , MS_O_k , MS_w , MS_O_m_bug + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - character(len) :: model_set (MS_model_count,2) = ' model ' , & - model_set_path = 'model_set.dat' , model_name , MS_columns(3)='' , & - MS_fields(30) , MS_parameters(30) + logical :: MS_recalculating = .false. - real(8) :: LumDist(N_models) = 1d0, MetrDist(N_models)= 1d0, DistMod(N_models)= 1d0, & !=- z for FCM > 0.0025 - dz, x_max, I1, RDR(2,N_RDR_grid), & + integer,parameter :: N_models = 9 , & + N_MS = 10 , & + MS_grid = 2d2 , & + MS_model_count = 1d2 - MS_table ( 1+3*MS_model_count , MS_grid ) = 0d0 + character(len) :: model_set (MS_model_count,2) = 'NaN' , model_name , MS_columns(3) , & + MS_fields(30) , MS_table_name = 'models.dat', MS_model(N_MS) = 'NaN' , & + MS_parameters(30) + + real(8) :: MS_z_max = 2d1 , MS_H_0 , MS_O_v , MS_O_k , MS_w , MS_O_m_bug , & + MS_Omega_w (N_MS) = NaN , MS_Omega_k(N_MS) = NaN , MS_EoS_w(N_MS) = NaN , & + MS_Hubble_0(N_MS) = NaN , MS_table ( 1+3*MS_model_count , MS_grid ) = 0d0 , & + + LumDist(N_models) = 1d0, MetrDist(N_models)= 1d0, DistMod(N_models)= 1d0 !=- z for FCM > 0.0025 + + integer :: MS_real_model_count = 0 , N_grid = 2d2 , & + MS_count_models = 0 , MS_count_w = 0 , MS_count_O_w = 0 , MS_count_O_k = 0 , MS_count_H_0 = 0 + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 contains + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine make_model_set + real(8) Omega_m + integer i,j,k,ii,ff + + if (model_set(1,1)=='NaN') then + call cheking_MS_parameters + + ii=0 + do jj=1,MS_count_models + + !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! + + if (MS_model(jj)=='wCDM') then + + do i= 1,MS_count_w + do j= 1,MS_count_O_w + do k= 1,MS_count_O_k + do kk=1,MS_count_H_0 + ii=1+ii + + if (MS_Omega_k(k)>0) then + ff=5 + else + ff=5 + endif + + Omega_m = abs ( 1 + MS_Omega_k(k) - MS_Omega_w(j) ) !=- the Friedmann equation + + !=- write(*,*) ' the model_set format:' + !=- write(*,*) ' wCDM H_0 w Omega_v Omega_k' + !=- write(*,*) ' FF H_0' + !=- write(*,*) ' TL H_0' + + model_set ( ii , 1 ) = & + ' wCDM ' // trim( realtostr( MS_Hubble_0(kk) ) ) // ' ' & + // trim( realtostr( MS_EoS_w(i) ) ) // ' ' & + // trim( realtostr( MS_Omega_w(j) ) ) // ' ' & + // trim( realtostr( MS_Omega_k(k) ) ) + + model_set ( ii , 2 ) = & + ' wCDM (w = ' // trim( realtostrf( MS_EoS_w(i) ,4 ,1) ) // & + ', {/Symbol W}_{DE} =' // trim( realtostrf( MS_Omega_w(j) ,4 ,1) ) // & + ', {/Symbol W}_m =' // trim( realtostrf( Omega_m ,4 ,1) ) // & + ', {/Symbol W}_k =' // trim( realtostrf( MS_Omega_k(k) ,ff,1) ) // & + ', H_0 =' // trim( realtostrf( MS_Hubble_0(kk),5 ,1) ) // ') ' - !=- write(*,*) ' the model_set format:' - !=- write(*,*) ' wCDM H_0 w Omega_v Omega_k' - !=- write(*,*) ' FF H_0' - !=- write(*,*) ' TL H_0' + enddo + enddo + enddo + enddo + + !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! + + else + + !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! - model_set ( 1 , 1 ) = ' wCDM 70 -1 0.7 0.0 ' - model_set ( 1 , 2 ) = ' the standard LCDM model ' + do kk=1,MS_count_H_0 + ii=1+ii + model_set ( ii , 1 ) = & + trim(MS_model(jj)) // ' ' // trim( realtostr( MS_Hubble_0(kk) ) ) - model_set ( 2 , 1 ) = ' wCDM 38.4 -1 0.15 -0.15 ' - model_set ( 2 , 2 ) = ' the dVMS model ' + model_set ( ii , 2 ) = trim(MS_model(jj)) & + // ', H_0 =' // trim( realtostrf( MS_Hubble_0(kk),5 ,1) ) // ') ' + end do - model_set ( 3 , 1 ) = ' wCDM 70 -1 1.0 0.0 ' - model_set ( 3 , 2 ) = ' the PV model ' + !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! - !model_set ( 7 , 1 ) = ' TL 70 ' - !model_set ( 7 , 2 ) = ' the TL model ' + endif + MS_real_model_count = ii + enddo - !model_set ( 6 , 1 ) = ' wCDM 70 -1 0.8 0.0 ' - !model_set ( 6 , 2 ) = ' the LCDM {/Symbol W}_{/Symbol L} = 0.8 model ' + else - model_set ( 4 , 1 ) = ' wCDM 70 -1 0.9 0.0 ' - model_set ( 4 , 2 ) = ' the LCDM {/Symbol W}_{/Symbol L} = 0.9 model ' + MS_real_model_count = count_noNaN_text(model_set(:,1)) - model_set ( 5 , 1 ) = ' FF 70 -1 0.9 0.0 ' - model_set ( 5 , 2 ) = ' the FF model ' + endif call compute_model_set call write_model_set call plot_model_set - end subroutine make_model_set + call MS_default + + end + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine cheking_MS_parameters + + MS_count_models = count_noNaN_text(MS_model) + MS_count_w = count_noNaN(MS_EoS_w) + MS_count_O_w = count_noNaN(MS_Omega_w) + MS_count_O_k = count_noNaN(MS_Omega_k) + MS_count_H_0 = count_noNaN(MS_Hubble_0) + + if ( MS_count_models==0 ) then + write(*,*) ' info: cheking_MS_parameters: MS_model = NaN and has been changed to LCDM' + MS_model(1) = 'wCDM' ; MS_count_models = 1 + endif + + if ( MS_count_w==0 ) then + write(*,*) ' info: cheking_MS_parameters: MS_EoS_w = NaN and has been changed to -1' + MS_EoS_w(1) = -1 ; MS_count_w = 1 + endif + + if ( MS_count_O_w==0 ) then + write(*,*) ' info: cheking_MS_parameters: MS_Omega_w = NaN and has been changed to 0.7' + MS_Omega_w(1) = O_7 ; MS_count_O_w = 1 + endif + + if ( MS_count_O_k==0 ) then + write(*,*) ' info: cheking_MS_parameters: MS_Omega_k = NaN and has been changed to 0.0' + MS_Omega_k(1) = 0d0 ; MS_count_O_k = 1 + endif + + if ( MS_count_H_0==0 ) then + write(*,*) ' info: cheking_MS_parameters: MS_Hubble_0 = NaN and has been changed to 70 km/s/Mpc' + MS_Hubble_0(1) = H70 ; MS_count_H_0 = 1 + endif + + end + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine MS_default + !=- for cosmology module + model_set (:,:) = 'NaN' + MS_table_name = 'models.dat' + MS_model(:) = 'NaN' + + MS_z_max = 2d1 + MS_Omega_w (:) = NaN ; MS_Omega_k(:) = NaN ; MS_EoS_w(:) = NaN + MS_Hubble_0(:) = NaN ; MS_table ( : , : ) = 0d0 + + LumDist(:) = 1d0 ; MetrDist(:)= 1d0 ; DistMod(:)= 1d0 !=- z for FCM > 0.0025 + + MS_real_model_count = 0 ; N_grid = 2d2 + MS_count_models = 0 ; MS_count_w = 0 ; MS_count_O_w = 0 ; MS_count_O_k = 0 ; MS_count_H_0 = 0 + + end subroutine MS_default + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine compute_model_set integer i,j ; MS_table(:,:)=0d0 - inquire( file = model_set_path , exist = file_exists ) + inquire( file = MS_table_name , exist = file_exists ) - if (.not. file_exists ) then + if (.not. file_exists .or. MS_recalculating ) then do i=1,MS_model_count @@ -144,7 +273,8 @@ subroutine compute_model_set dm_TL( MS_table ( 1 , j ) , MS_H_0 ) case default - write(*,*) 'MS_compute: unknown model mask: ' // trim(model_name) + if (model_name/='NaN') & + write(*,*) 'MS_compute: unknown model mask: ' // trim(model_name) end select enddo @@ -158,21 +288,21 @@ subroutine compute_model_set endif end subroutine compute_model_set - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine write_model_set integer i character(len) title(MS_real_model_count) - inquire( file = model_set_path , exist = file_exists ) + inquire( file = MS_table_name , exist = file_exists ) - if (MS_model_count>0 .and. .not. file_exists ) then + if (MS_model_count>0 .and. .not. file_exists .or. MS_recalculating ) then MS_columns(1) = ':metric_r' MS_columns(2) = ':luminosity_d_L' MS_columns(3) = ':distance_modulus' - open(unit_1,file=model_set_path,status='replace') + open(unit_1,file=MS_table_name,status='replace') line=' ' do i=1,MS_real_model_count @@ -198,7 +328,7 @@ subroutine write_model_set endif end subroutine write_model_set - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! integer function fix_index( in_index , max_border ) integer in_index , out_index , max_border @@ -209,30 +339,35 @@ integer function fix_index( in_index , max_border ) fix_index = out_index end function fix_index - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine plot_model_set + character(len) logscale_interlineation , log_y , logscale_lable integer i,j,k call clear_plot_fields + pngterm = 'set term pngcairo enhanced font "Verdana,20" size 1400, 1050' + epsterm = 'set term postscript enhanced color font "Verdana,14" size 8.5, 6.3' + GNUfields(logscale) = 'set logscale' - GNUfields(format_y) = 'set format y "10^{%L}"' - GNUfields(xrange) = 'set xrange [0.01:*]' + GNUfields(format_y) = '10^{%L}' + GNUfields(xrange) = '[0.01:*]' + + GNUfields(legend) = 'right bottom' - GNUfields(legend) = 'set key right bottom' - GNUfields(xlabel) = 'set xlabel "redshift z"' + GNUfields(grid) = 'xtics ytics mxtics mytics' MS_fields(1) = 'r(z), [Mpc]' MS_fields(2) = 'd_L(z), [Mpc]' MS_fields(3) = '{/Symbol m}(z)' - MS_fields(4) = 'log10 r(z)/r_{model-1}(z)' - MS_fields(5) = 'log10 d_L(z)/d_{L,model-1}(z)' - MS_fields(6) = 'log10 {/Symbol m}(z) - mu_{model-1}(z)' + MS_fields(4) = 'log10 r(z)/r_{' // trim( model_set (1,2) ) // '}(z)' + MS_fields(5) = 'log10 d_L(z)/d_{L,' // trim( model_set (1,2) ) // '}(z)' + MS_fields(6) = '{/Symbol m}(z) - {/Symbol m}_{' // trim( model_set (1,2) ) // '}(z)' MS_fields(7) = 'log10 {/Symbol D}r(z)' MS_fields(8) = 'log10 {/Symbol D}d_L(z)' - MS_fields(9) = 'log10 {/Symbol D}{/Symbol m}(z)' + MS_fields(9) = '{/Symbol D}{/Symbol m}(z)' MS_fields(10) = 'r(z)' MS_fields(11) = 'd_L(z)' @@ -242,48 +377,72 @@ subroutine plot_model_set MS_fields(14) = 'delta_d_L(z)' MS_fields(15) = 'delta_mu(z)' - do k=1,2 - do j=1,3 + do k=1,2 !=- Delta + do j=1,3 !=- r(z), d_L(z), mu(z) - GNUfields(title) = 'set title "'// trim(MS_fields(j+6*(k-1))) //' for different models"' + GNUfields(mxtics) = '2' + GNUfields(mytics) = '2' + logscale_lable='' + if ( k==1 ) then + logscale_lable = ', logscale' + GNUfields(mxtics) = '5' + GNUfields(mytics) = '5' + endif + GNUfields(xlabel) = 'set xlabel "redshift z' // trim(logscale_lable) // '"' + if (j==3) logscale_lable = '' + + GNUfields(title) = 'set title "'// trim(MS_fields(j+6*(k-1))) //' for different models (CosMod v1.0)"' + + logscale_interlineation='-logscale' + log_y='log10' + if (k==2) logscale_interlineation = '' + if (k==2) GNUfields(legend) = 'set key off' if (j==3 .and. k/=2) then - GNUfields(logscale) = '#set logscale' + GNUfields(logscale) = 'set logscale x' GNUfields(format_y) = '#set format y "10^{%L}"' endif + if (k==2 .and. j==3) log_y='' + + + if (MS_real_model_count==0) MS_real_model_count = count_noNaN_text(model_set (:,1)) do i=1,MS_real_model_count select case (k) case (1) - GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(j+6*(k-1))) // '"' + GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(j+6*(k-1))) // trim(logscale_lable) // '"' + + GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(MS_table_name)) // & + '" u 1:' // trim(inttostr(1+j+3*(i-1))) // ' w l ls ' // trim(inttostr(i)) // ' lw 3 ' - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:' // trim(inttostr(1+j+3*(i-1))) // ' w l ls ' // trim(inttostr(i)) GNUfields (title1+(i-1)*2) = ' title "' // trim( model_set (i,2) ) // '"' case(2) - GNUfields(yrange) = 'set yrange [-0.5:1]' - if (j==3) GNUfields(yrange) = 'set yrange [-0.025:0.05]' + GNUfields(yrange) = 'set yrange [*:*]' + !=- GNUfields(yrange) = 'set yrange [-0.5:1]' + !=- if (j==3) GNUfields(yrange) = 'set yrange [-0.025:0.05]' GNUfields(ylabel) = 'set ylabel "' // trim(MS_fields(j+6*(k-1))) // ' = ' // & trim(MS_fields(j+3*(k-1))) // '"' - GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(model_set_path)) // & - '" u 1:(log10($' // trim(inttostr(1+j+3*(i-1))) // '/$' // & - trim(inttostr(1+j)) // ')) w l ls ' // trim(inttostr(i)) + GNUfields (plot1+(i-1)*2) = ', "' // trim(slashfix(MS_table_name)) // & + '" u 1:('// trim(log_y) //'($' // trim(inttostr(1+j+3*(i-1))) // ') - '// & + trim(log_y) //'($' // & + trim(inttostr(1+j)) // ')) w l ls ' // trim(inttostr(i)) // ' lw 3 ' + GNUfields (title1+(i-1)*2) = ' title "' // trim( model_set (i,2) ) // '"' end select end do GNUfields (plot1) = GNUfields (plot1)(2:len) - !graph_name = MS_fields(j+9+3*(k-1)) // '-logscale' - ! GNUfields(extention_out_figure)='#eps' ; call plot(model_set_path) - graph_name = MS_fields(j+9+3*(k-1)) // '-logscale' - GNUfields(extention_out_figure)='#png' ; call plot(model_set_path) + graph_name = trim(MS_fields(j+9+3*(k-1))) // trim(logscale_interlineation) + GNUfields(extention_out_figure)='#eps' ; call plot(MS_table_name) + graph_name = trim(MS_fields(j+9+3*(k-1))) // trim(logscale_interlineation) + GNUfields(extention_out_figure)='#png' ; call plot(MS_table_name) enddo enddo end subroutine plot_model_set - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine redshift_distance_relation !=- RDR( redshift , luminosity distance ) if (RDR_calculating) then @@ -296,6 +455,8 @@ subroutine redshift_distance_relation !=- RDR( redshift , luminosity distance ) endif end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function fun_from_R_to_z(luminosity_distance) real(8) luminosity_distance ; fun_from_R_to_z = 0d0 call redshift_distance_relation @@ -304,6 +465,8 @@ real(8) function fun_from_R_to_z(luminosity_distance) end do end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function fun_from_z_to_R(redshift,H0) real(8) redshift,H0 ; fun_from_z_to_R = 0d0 select case (RDR_type) @@ -318,28 +481,36 @@ real(8) function fun_from_z_to_R(redshift,H0) end select end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function mag_to_Lum( visible_magnitude , luminosity_distance ) real(8) visible_magnitude,luminosity_distance ; Mag_to_Lum = 0d0 if (visible_magnitude .ne. 0) Mag_to_Lum = L_sun * 1d1**( 0.4d0*( M_sun - & vis_to_abs_mag( visible_magnitude , luminosity_distance ) ) ) end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function Lum_to_Mag(Luminosity) real(8) Luminosity ; Lum_to_Mag = 0d0 Lum_to_Mag = M_sun - 2.5d0 * log10 ( Luminosity / L_sun ) end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function abs_to_vis_mag( absolute_magnitude , luminosity_distance ) real(8) absolute_magnitude,luminosity_distance ; abs_to_vis_mag = 0d0 if (absolute_magnitude .ne. 0) abs_to_vis_mag = absolute_magnitude - 5d0 + 5d0 * log10 ( luminosity_distance * 1d5 ) !=- 10 pc end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function vis_to_abs_mag( visible_magnitude , luminosity_distance ) real(8) visible_magnitude,luminosity_distance ; vis_to_abs_mag = 0d0 if (visible_magnitude .ne. 0) vis_to_abs_mag = visible_magnitude + 5d0 - 5d0 * log10 ( luminosity_distance * 1d5 ) !=- 10 pc end function - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! !real(8) function R_LCDM(z) ; real(8) z ; R_LCDM = dl * IntHwCDM(z,-1d0,O_m,0d0) ; end function @@ -373,6 +544,8 @@ real(8) function dm_FCM(z,H0) ; real(8) z,H0 ; dm_FCM = mu( D_FCM(z,H0) real(8) function dm_MM(z,H0) ; real(8) z,H0 ; dm_MM = mu( D_MM(z,H0) ) ; end function !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=! + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function YW(z) real(8) z YW=0d0 ; N=1d4 ; x_max=6d0 @@ -406,6 +579,8 @@ real(8) function YW(z) endif end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function IntHwCDM(z,w,Om,Ok) integer i,N real(8) z,w,Om,Ok,Ow @@ -436,237 +611,8 @@ real(8) function RzwCDM(z,w,Om,Ok) end if end if end function RzwCDM -!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 - subroutine models_head(titles) - character(len) titles(N_titles) - titles(1) ='z ' ; - titles(2) ='z+1' ; - titles(3) ='log(z+1)' ; - - titles(4) ='dLCDM/rLCDM' ; titles(22) ='dLCDM' - titles(5) ='d_CSS(z)/rLCDM' ; titles(23) ='d_CSS(z)' - titles(6) ='d_TL(z)/rLCDM' ; titles(24) ='d_TL(z)' - titles(7) ='d_FCM(z)/rLCDM' ; titles(25) ='d_FCM(z)' - titles(8) ='d_MM(z)/rLCDM' ; titles(26) ='d_MM(z)' - titles(9) ='d_wCDM(z,-1.0,0.3,-0.1)/rLCDM' ; titles(27) ='d_wCDM(z,-1.0,0.3,-0.1)' - titles(10) ='d_wCDM(z,-1.0,0.3, 0.1)/rLCDM' ; titles(28) ='d_wCDM(z,-1.0,0.3, 0.1)' - titles(11) ='d_wCDM(z,-1.0,0.1, 0.0)/rLCDM' ; titles(29) ='d_wCDM(z,-1.0,0.1, 0.0)' - titles(12) ='d_wCDM(z,-0.5,0.5, 0.2)/rLCDM' ; titles(30) ='d_wCDM(z,-0.5,0.5, 0.2)' - - titles(13) ='rLCDM /rLCDM' ; titles(31) ='rLCDM' - titles(14) ='r_CSS(z)/rLCDM' ; titles(32) ='r_CSS(z)' - titles(15) ='r_TL(z)/rLCDM' ; titles(33) ='r_TL(z)' - titles(16) ='r_FCM(z)/rLCDM' ; titles(34) ='r_FCM(z)' - titles(17) ='r_MM(z)/rLCDM' ; titles(35) ='r_MM(z)' - titles(18) ='r_wCDM(z,-1.0,0.3,-0.1)/rLCDM' ; titles(36) ='r_wCDM(z,-1.0,0.3,-0.1)' - titles(19) ='r_wCDM(z,-1.0,0.3, 0.1)/rLCDM' ; titles(37) ='r_wCDM(z,-1.0,0.3, 0.1)' - titles(20) ='r_wCDM(z,-1.0,0.1, 0.0)/rLCDM' ; titles(38) ='r_wCDM(z,-1.0,0.1, 0.0)' - titles(21) ='r_wCDM(z,-0.5,0.5, 0.2)/rLCDM' ; titles(39) ='r_wCDM(z,-0.5,0.5, 0.2)' - - titles(40) ='dm_LCDM-dm_LCDM' ; titles(49) ='dm_LCDM' - titles(41) ='dm_CSS(z)-dm_LCDM' ; titles(50) ='dm_CSS(z)' - titles(42) ='dm_TL(z)-dm_LCDM' ; titles(51) ='dm_TL(z)' - titles(43) ='dm_FCM(z)-dm_LCDM' ; titles(52) ='dm_FCM(z)' - titles(44) ='dm_MM(z)-dm_LCDM' ; titles(53) ='dm_MM(z)' - titles(45) ='dm_wCDM(z,-1.0,0.3,-0.1)-dm_LCDM' ; titles(54) ='dm_wCDM(z,-1.0,0.3,-0.1)' - titles(46) ='dm_wCDM(z,-1.0,0.3, 0.1)-dm_LCDM' ; titles(55) ='dm_wCDM(z,-1.0,0.3, 0.1)' - titles(47) ='dm_wCDM(z,-1.0,0.1, 0.0)-dm_LCDM' ; titles(56) ='dm_wCDM(z,-1.0,0.1, 0.0)' - titles(48) ='dm_wCDM(z,-0.5,0.5, 0.2)-dm_LCDM' ; titles(57) ='dm_wCDM(z,-0.5,0.5, 0.2)' - - end subroutine - - subroutine models_calculating(datafile) - character(len) datafile - - open(1,file=datafile, status='replace'); call models_head(titles) - theformat='(A2,'//trim(inttostr(N_titles))//'(i2,1x,A22,1x))' - do j=1,N_titles - write(1,'(A2,i2,1x,A40)') '# ',j,titles(j) - enddo - write(1,theformat) '# ',(j,titles(j),j=1,N_titles) - - do i=1,N_grid - z = dsqrt(z_max)/N_grid * i!(i+0.01d0) - z = z*z - - LumDist(1) = d_LCDM (z,H_0) ; MetrDist(1) = r_LCDM (z,H_0) - LumDist(2) = d_CSS (z,H_0) ; MetrDist(2) = r_CSS (z,H_0) - LumDist(3) = d_TL (z,H_0) ; MetrDist(3) = r_TL (z,H_0) - LumDist(4) = d_FCM (z,H_0) ; MetrDist(4) = r_FCM (z,H_0) - LumDist(5) = d_MM (z,H_0) ; MetrDist(5) = r_MM (z,H_0) - LumDist(6) = d_wCDM (z,H_0,-1.0d0,0.85d0, 0.0d0) ; MetrDist(6) = r_wCDM (z,H_0,-1.0d0,0.85d0, 0.0d0) - LumDist(7) = d_wCDM (z,H_0,-1.0d0,0.20d0, 0.0d0) ; MetrDist(7) = r_wCDM (z,H_0,-1.0d0,0.20d0, 0.0d0) - LumDist(8) = d_wCDM (z,H_0,-1.0d0,0.1d0, 0.0d0) ; MetrDist(8) = r_wCDM (z,H_0,-1.0d0,0.1d0, 0.0d0) - LumDist(9) = d_wCDM (z,H_0,-0.5d0,0.5d0, 0.2d0) ; MetrDist(9) = r_wCDM (z,H_0,-0.5d0,0.5d0, 0.2d0) - - DistMod(1) = dm_LCDM (z,H_0) - DistMod(2) = dm_CSS (z,H_0) - DistMod(3) = dm_TL (z,H_0) - DistMod(4) = dm_FCM (z,H_0) - DistMod(5) = dm_MM (z,H_0) - DistMod(6) = dm_wCDM (z,H_0,-1.0d0,0.85d0, 0.0d0) - DistMod(7) = dm_wCDM (z,H_0,-1.0d0,0.20d0, 0.0d0) - DistMod(8) = dm_wCDM (z,H_0,-1.0d0,0.1d0, 0.0d0) - DistMod(9) = dm_wCDM (z,H_0,-0.5d0,0.5d0, 0.2d0) - - do j=1,N_models; if ( LumDist(j)==0) LumDist(j)=1d0 - if (MetrDist(j)==0) MetrDist(j)=1d0 ; end do - - theformat2='('//trim(inttostr(N_titles))//'(F20.8,6x))' - write(1,theformat2) z,z+1,log(z+1) , & - ( LumDist(j)/ LumDist(1),j=1,N_models) , & - (MetrDist(j)/MetrDist(1),j=1,N_models) , & - ( LumDist(j) ,j=1,N_models) , & - (MetrDist(j) ,j=1,N_models) , & - ( DistMod(j)- DistMod(1),j=1,N_models) , & - ( DistMod(j) ,j=1,N_models) - end do - - close(1) - - end subroutine - -!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 - subroutine models(datafile) - character(len) datafile - - call models_calculating(datafile) - - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - GNUfields(1) = 'set linestyle 1 lw 3 pt 7 ps 0.9 lt rgb "red"' - GNUfields(2) = 'set linestyle 2 lw 3 pt 7 ps 0.9 lt rgb "grey"' - GNUfields(3) = 'set linestyle 3 lw 3 pt 7 ps 0.9 lt rgb "black"' - GNUfields(4) = 'set linestyle 4 lw 3 pt 7 ps 0.9 lt rgb "green"' - - GNUfields(5) = 'set linestyle 5 lw 3 pt 7 ps 0.9 lt rgb "brown"' - GNUfields(22) = 'set linestyle 6 lw 3 pt 7 ps 0.9 lt rgb "blue"' - GNUfields(23) = 'set linestyle 9 lw 3 pt 7 ps 0.9 lt rgb "brown"' - - GNUfields(10) = 'set yrange[-0.6:1.0]' - GNUfields(9) = 'set xlabel "z"' ; j = 1 - GNUfields(14) = 'set ylabel "{/Symbol D}log d_L(z) = log d_L(z)/d_L_{,{/Symbol L}-CDM}(z)"' - GNUfields(7) = 'set title "{/Symbol D}log d_L(z) for different models"' - !GNUfields(24) = 'set key opaque' - - !GNUfields(31) = 'set key left bottom' - GNUfields(31) = 'set nokey' - - do i=1,N_models - GNUfields(31+i+i) = ', "" u '//trim(inttostr(j))//':(log10($'//trim(inttostr(i+3+N_models*0))// & - ')) w l ls '//trim(inttostr(i)) - enddo - - GNUfields(34) = ' title "{{/Symbol L}-CDM}({/Symbol W}_{/Symbol L}=0.7)"' - GNUfields(36) = ' title "{{/Symbol L}-CDM}({/Symbol W}_{/Symbol L}=1.0),CSS"' - GNUfields(38) = ' title "{TL}"' - GNUfields(40) = ' title "{FF}"' - GNUfields(42) = ' title "{{/Symbol L}-CDM}({/Symbol W}_{/Symbol L}=0.9)"' - GNUfields(44) = ' title " w-CDM(w=-0.5,{/Symbol W}_{/Symbol L}=0.7,{/Symbol W}_k=0.2)"' - - GNUfields(20) = '#set logscale y' - GNUfields(12) = '#' - GNUfields(17) = '#' - GNUfields(11) = 'set mxtics 2' - GNUfields(16) = 'set mytics 2' - - graph_name = 'Shirokov_fig6' ; GNUfields(extention_out_figure)='#eps' ; call plot(datafile) - graph_name = '_d_L(z)' ; GNUfields(extention_out_figure)='#png' ; call plot(datafile) - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - GNUfields(14) = 'set ylabel "{/Symbol D}log r(z) = log r(z)/r_{{/Symbol L}-CDM}(z)"' - GNUfields(7) = 'set title "{/Symbol D}log r(z) for different models"' - - do i=1,N_models - GNUfields(31+i+i) = ', "" u '//trim(inttostr(j))//':(log10($'//trim(inttostr(i+3+N_models*1))// & - ')) w l ls '//trim(inttostr(i)) - enddo - - graph_name = 'Shirokov_fig4' ; GNUfields(extention_out_figure)='#eps' ; call plot(datafile) - graph_name = '_d_r(z)' ; GNUfields(extention_out_figure)='#png' ; call plot(datafile) - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - GNUfields(11) = 'set mxtics 5' - GNUfields(16) = 'set mytics 5' - - GNUfields(9) = 'set xlabel "z"' ; ii = 1 - !GNUfields(10) = 'set xrange[0.1:2]' - !GNUfields(15) = 'set yrange[100:10000]' - GNUfields(10) = 'set xrange[1:20]' - GNUfields(15) = 'set yrange[500:1000000]' - GNUfields(20) = 'set logscale' - - GNUfields(17) = 'set format y "10^{%L}"' - GNUfields(14) = 'set ylabel "d_L(z), [Mpc]"' - GNUfields(7) = 'set title "d_L(z) for different models in logscales"' - - GNUfields(31) = 'set key right bottom' - - do i=1,N_models - GNUfields(31+i+i) = ', "" u '//trim(inttostr(j))//':(($'//trim(inttostr(i+3+N_models*2))// & - ')) w l ls '//trim(inttostr(i)) - enddo - - graph_name = 'Shirokov_fig5' ; GNUfields(extention_out_figure)='#eps' ; call plot(datafile) - graph_name = '_d_L(z)-logscale' ; GNUfields(extention_out_figure)='#png' ; call plot(datafile) - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - GNUfields(14) = 'set ylabel "r(z), [Mpc]"' - GNUfields(7) = 'set title "r(z) for different models in logscales"' - GNUfields(15) = 'set yrange[500:100000]' - - do i=1,N_models - GNUfields(31+i+i) = ', "" u '//trim(inttostr(j))//':(($'//trim(inttostr(i+3+N_models*3))// & - ')) w l ls '//trim(inttostr(i)) - enddo - - graph_name = 'Shirokov_fig3' ; GNUfields(extention_out_figure)='#eps' ; call plot(datafile) - graph_name = '_r(z)-logscale' ; GNUfields(extention_out_figure)='#png' ; call plot(datafile) - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - GNUfields(14) = 'set ylabel "{/Symbol D}{/Symbol m}(z) = {/Symbol m}(z) - {/Symbol m}_{{/Symbol L}-CDM}(z)"' - GNUfields(7) = 'set title "{/Symbol D}{/Symbol m}(z) for different models"' - GNUfields(15) = 'set yrange[-3:5]' - GNUfields(10) = 'set xrange[0.001:20]' - GNUfields(20) = 'unset logscale' - GNUfields(17) = 'set format y "%g"' - GNUfields(11) = 'set mxtics 2' - GNUfields(16) = 'set mytics 2' - - GNUfields(31) = 'set nokey' - - do i=1,N_models - GNUfields(31+i+i) = ', "" u '//trim(inttostr(j))//':(($'//trim(inttostr(i+3+N_models*4))// & - ')) w l ls '//trim(inttostr(i)) - enddo - - GNUfields(34) = ' title "{/Symbol m}_{{/Symbol L}-CDM}({/Symbol W}_{/Symbol L}=0.7)"' - GNUfields(36) = ' title "{/Symbol m}_{{/Symbol L}-CDM}({/Symbol W}_{/Symbol L}=1.0),CSS"' - GNUfields(38) = ' title "{/Symbol m}_{TL}"' - GNUfields(40) = ' title "{/Symbol m}_{FF}"' - GNUfields(42) = ' title "{/Symbol m}_{{/Symbol L}-CDM}({/Symbol W}_{/Symbol L}=0.9)"' - GNUfields(44) = ' title " {/Symbol m}_{w-CDM}(w=-0.5,{/Symbol W}_{/Symbol L}=0.7,{/Symbol W}_k=0.2)"' - - graph_name = 'Shirokov_fig8' ; GNUfields(extention_out_figure)='#eps' ; call plot(datafile) - graph_name = '_u(z)' ; GNUfields(extention_out_figure)='#png' ; call plot(datafile) - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - GNUfields(14) = 'set ylabel "{/Symbol m}(z)"' - GNUfields(7) = 'set title "{/Symbol m}(z) for different models in z-logscale"' - GNUfields(15) = 'set yrange[32:57]' - GNUfields(10) = 'set xrange[0.1:20]' - GNUfields(20) = 'set logscale x' - GNUfields(11) = 'set mxtics 5' - GNUfields(16) = 'set mytics 5' - - GNUfields(31) = 'set key right bottom' - - do i=1,N_models - GNUfields(31+i+i) = ', "" u '//trim(inttostr(j))//':(($'//trim(inttostr(i+3+N_models*5))// & - ')) w l ls '//trim(inttostr(i)) - enddo - - graph_name = 'Shirokov_fig7' ; GNUfields(extention_out_figure)='#eps' ; call plot(datafile) - graph_name = '_u(z)-logscale' ; GNUfields(extention_out_figure)='#png' ; call plot(datafile) - !-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - - - end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! end module !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/GNUplot.f95 b/GNUplot.f95 index 5167fc7..83e8eac 100644 --- a/GNUplot.f95 +++ b/GNUplot.f95 @@ -7,14 +7,16 @@ module GNUplot !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 integer, parameter :: N_GNUfields = 200 - integer :: ls1=1, ls2=2, ls3=3, ls4=4, ls5=5, title=7,xlabel=9,xrange=10,mxtics=11,format_x=12, & + integer :: ls1=1, ls2=2, ls3=3, ls4=4, ls5=5, title=6,xlabel=9,xrange=10,mxtics=11,format_x=12, & ylabel=14,yrange=15,mytics=16,format_y=17,grid=19,logscale=20,add1=23,add2=24,add3=25, add4=26, add5=27, add6=28, & extention_out_figure=29,legend=31,plot1=33,title1=34,plot2=35,title2=36,plot3=37,title3=38,plot4=39,title4=40, & plot5=41,title5=42,plot6=43,title6=44,plot7=45,title7=46,plot8=47,title8=48,plot9=49,title9=50,plot10=51,title10=52, & plot11=53,title11=54,plot12=55,title12=56,plot13=57,title13=58,plot14=59,title14=60,set_parametric=22, & thickness , comma_fix , dash_type - character(len) :: pngterm = 'set term pngcairo enhanced font "Verdana,10" size 850, 630' , & + real(8) :: external_parameter = 0 + + character(len) :: pngterm = 'set term pngcairo enhanced font "Verdana,12" size 1400, 1050' , & epsterm = 'set term postscript enhanced color font "Verdana,14" size 8.5, 6.3' , & scriptpath = '', figurepath = '', figoutput = '', GNUdatafile = '', graph_name='', & fig_dir = '', GNUplots = '', GNUscripts = '', script_file_path='', plot_dir='', & @@ -34,7 +36,7 @@ module GNUplot ! ' set ylabel "y" # 14 y-axis label ', & ! ' set yrange [0:1] # 15 ', & ! ' set mytics 5 # 16 ', & -! ' set format y "10^{%L}" # 17 ', & +! ' set format y "%2.0t{/Symbol \327}10^{%L}" # 17 ', & ! ' ', & ! ' set grid xtics ytics mxtics mytics # 19 ', & ! ' set logscale # 20 ', & @@ -64,28 +66,24 @@ subroutine plot(datafile) !-=- plot autoscripting unit_2 = random_unit() - select case (operation_system) + select case (operating_system) case(1) !=- Windows datafile = backslashfix(datafile) GNUplots = 'plots\' GNUscripts = 'plots\GNUscripts\' plot_dir = trim(file_path(datafile)) // trim(GNUplots) // trim(name(datafile)) & // '\' // trim(fig_dir) - system_commands(1) = 'MD' - system_commands(2) = 'move' case(2) !=- Linux datafile = slashfix(datafile) GNUplots = 'plots/' GNUscripts = 'plots/GNUscripts/' plot_dir = trim(file_path(datafile)) // trim(GNUplots) // trim(name(datafile)) & // '/' // trim(fig_dir) - system_commands(1) = 'mkdir -p' - system_commands(2) = 'mv' end select - call system( trim(system_commands(1)) // ' ' // trim(file_path(datafile)) // & - trim(GNUscripts) // ' >> log.log ' ) - call system( trim(system_commands(1)) // ' ' // trim(plot_dir) // ' >> log.log ' ) + output = trim(file_path(datafile)) // trim(GNUscripts) + call create_folder( output ) + call create_folder( plot_dir ) figoutput = trim(file_path(datafile)) // trim(GNUscripts) // trim(name(datafile)) // & '_' // trim(graph_name) @@ -98,7 +96,7 @@ subroutine plot(datafile) !-=- plot autoscripting plot_extension = '.eps' case default write(unit_2,*) pngterm - plot_extension = '.png' + plot_extension = '.png' end select figoutput = trim(figoutput) // trim(plot_extension) @@ -109,12 +107,13 @@ subroutine plot(datafile) !-=- plot autoscripting write(unit_2,'(A)') trim(GNUfields(i)) enddo !write(9,*) 'plot "> log.log ' ) + call move_file( figoutput , plot_dir ) graph_name='' @@ -136,6 +135,13 @@ end function GNUcorrect subroutine GNUfields_fix + character(len) plot + + read(GNUfields(plot1),*) plot + if ( plot/='plot' .and. plot/='splot' ) GNUfields(plot1) = 'plot ' // trim(GNUfields(plot1)) + + if (GNUfields(extention_out_figure)(1:1)/='#') & + GNUfields(extention_out_figure) = '#' // trim(GNUfields(extention_out_figure)) if ( .not. GNUcorrect( GNUfields( title ) ) ) & GNUfields( title ) = 'set title "' // trim(GNUfields( title )) // '"' @@ -167,7 +173,7 @@ subroutine GNUfields_fix if ( GNUfields(i) /= '' .and. .not. word_search(GNUfields(i),'title') & .and. .not. word_search(GNUfields(i),'notitle') ) & - GNUfields(i) = ' title "' // trim(GNUfields(i)) // '"' + GNUfields(i) = ' title "' // trim(GNUfields(i)) // '"' end do if ( .not. GNUcorrect( GNUfields( format_y ) ) ) & @@ -184,9 +190,11 @@ subroutine clear_plot_fields GNUfields( : ) = '' !=- default fig_dir = '' - pngterm = 'set term pngcairo enhanced font "Verdana,10" size 850, 630' !=- default + pngterm = 'set term pngcairo enhanced font "Verdana,16" size 1400, 1050' epsterm = 'set term postscript enhanced color font "Verdana,14" size 8.5, 6.3' + !=- Verdana , Palatino-Roman + end subroutine diff --git a/global.f95 b/global.f95 index c61e625..b33140d 100644 --- a/global.f95 +++ b/global.f95 @@ -5,24 +5,27 @@ module global implicit none - integer :: operation_system = 0 !=- 1 Windows, 2 Linux + integer :: operating_system = 0 !=- 1 Windows, 2 Linux - logical :: file_exists + logical :: file_exists , index_info , write_percent_fix , all_error_reporst = .false. integer,parameter :: len=256, N_col_GRB = 90, N_GRB=7d3, & !=- constants - N_titles = 57, N_comparisons = 4, N_files = 20, N_folders = 10 !=- 3+6*N_models + N_titles = 57, N_comparisons = 4, N_files = 20, N_folders = 100, & !=- 3+6*N_models + N_commands = 10 integer :: i,ii,iii,j,jj,jjj,k,kk,kkk,n,nn,nnn,m,mm,mmm, & unit_1 = 1, unit_2 = 2, unit_3 = 3 , unit_4 = 4 , unit_5 = 5 , unit_6 = 6 , unit_7 = 7 , & - counter , final_counter + counter , final_counter , iostat_value , i_percent , N_percent , mod_percent = 10 , m_percent character(len) :: theformat='', theformat2='', datafile='', titles(N_titles)='', & input='', output='', line='',preformat='',path='',head_format='', & text1='',text2='',text3='',text4='',text5='', command='',figure_number='', & - files(N_files)='',new_files(N_files)='', folders(N_folders), & - system_commands(10) + files(N_files)='NaN',new_files(N_files)='NaN', folders(N_folders) = 'NaN', & + system_commands(10), filepath='', columns, filename='', commandN(N_commands) - real(8) a,b,bbb,d,e,g,l,t,tx,ty,nx,ny,aa,bb,vva,vvb,q,p,sa,sb,ssa,ssb,sm, & + real(8),parameter :: NaN = 1d150 + + real(8) a,b,bbb,d,e,g,l,t,tx,ty,nx,ny,aa,bb,vva,vvb,q,p,sa,sb,ssa,ssb,sm, lb,rb, & x,xx,xxx, y,yy,yyy, z,zz,zzz, w,ww,www, v,vv,vvv,r, maximum, & alpha,beta,p1,p2,p3,p4,p5,sumx,sumx2,sumy,sumxy , & z_max_var, log_z, GRB_shift , & @@ -38,20 +41,36 @@ subroutine define_system call system(' echo %OS% > system.log') line = read_last_string('system.log') - if ( line == 'Windows_NT' ) operation_system = 1 + if ( line == 'Windows_NT' ) operating_system = 1 call system(' echo $SHELL >> system.log') line = read_last_string('system.log') - if ( line == '/bin/bash' ) operation_system = 2 + if ( line == '/bin/bash' ) operating_system = 2 + select case (operating_system) + case(1) !=- Windows + write(*,*) 'OS: Windows_NT' + call sleep (1) + call system('cls') + case(2) !=- Linux + write(*,*) 'OS: Linux' + call sleep (1) + !call system('clear') + case default !=- Other + write(*,*) 'define_system: wrong mask' + call sleep (10) ; operating_system = 1 + end select end subroutine define_system integer function random_unit() real(8) x - call random_number(x) - random_unit = int(1d3*x) + do + call random_number(x) + random_unit = int(1d3*x) + if (random_unit>1d2) exit + enddo end function @@ -68,9 +87,30 @@ character(len) function read_last_string(log_pathfile) + subroutine prepare_percent( N_max ) + integer N_max + i_percent=0 ; N_percent = N_max ; call CPU_TIME(start_time) ; m_percent=0 + end subroutine prepare_percent + + + + subroutine write_percent !=- before a cycle set `call prepare_percent( N_max )` + i_percent = 1+i_percent + a = 1d2*i_percent/N_percent + if ( int(a)>m_percent .and. int(a)>0 .and. (mod(int(a),mod_percent) == 0 .or. int(a) == 1 .or. int(a) == 99) ) then + m_percent=int(a) + call CPU_TIME(final_time) + b = final_time - start_time + write(*,'(6x,"complited",i4,"%, work time is",F10.2,"s, time left is",F10.2,"s, processed points are",i10," from",i10)') & + int(a),b,b/int(a)*1d2-b, i_percent,N_percent + endif + end subroutine write_percent + + + character(len) function make_workdir() - select case (operation_system) + select case (operating_system) case(1) !=- Windows call system( 'echo %cd% > log.log' ) @@ -84,7 +124,7 @@ character(len) function make_workdir() make_workdir = file_path(line) case default !=- Other - write(*,*) 'fatal error: make_workdir: unknown mask' + write(*,*) 'fatal error: make_workdir: unknown mask' end select end function make_workdir @@ -92,20 +132,36 @@ end function make_workdir character(len) function name(path) - character(len) path; call file_name(path,ii,jj); name=''; name=path(ii:jj) + character(len) path; integer ii,jj + call file_name(path,ii,jj); name=''; name=path(ii:jj) + end function + + character(len) function name_ext(path) + character(len) path; integer ii,jj + call file_name(path,ii,jj); name_ext=''; name_ext=path(ii:len) end function subroutine file_name(path,l1,l2) !-=- file path excision - character(len) path ; integer i,l1,l2 + character(len) path ; integer i,l1,l2 ; l1=1 ; l2=len do i=len,1,-1 ; if (path(i:i)=='.') then ; l2=i-1 ; exit ; endif ; enddo do i=len,1,-1 ; if (path(i:i)=='\'.or.path(i:i)=='/') then ; l1=i+1 ; exit ; endif ; enddo end subroutine file_name character(len) function file_path(path) !-=- file name excision - integer i - character(len) path ; integer l1,l2 - do i=len,1,-1 ; if (path(i:i)=='\'.or.path(i:i)=='/') then ; l1=i ; exit ; endif ; enddo - file_path='' ; file_path = path(1:l1) + character(len) path ; integer l1,l2,i + file_path = path ; l2=len + do i=len,1,-1 + if (path(i:i)=='\'.or.path(i:i)=='/') then + l2=i + file_path = path(1:l2) + !write(*,*) l1, trim(file_path) + exit + endif + enddo + !call file_name(path,l1,l2) + !write(*,*) l1 , l2 , (l2==len) , trim(path) , ' ' , trim(name(path)) + if (l2==len) file_path = '' ! trim(name(path)) + !write(*,*) l1 , l2 , (l2==len) , trim(file_path) , ' ' , trim(name(file_path)) end function character(len) function inttostr(integer_number) !=- Convert an integer to string @@ -137,8 +193,10 @@ character(len) function realtostrff(real_number,theformat) !=- INVALID end function realtostrff integer function file_volume(file_path) - character(len) file_path ; file_volume=0 - unit_4 = random_unit() + character(len) file_path + integer unit_4 + + unit_4 = random_unit() ; file_volume=0 inquire( file = file_path , exist = file_exists ) if ( file_exists ) then @@ -150,7 +208,7 @@ integer function file_volume(file_path) 11 continue ; close(unit_4) else write(*,*) ' file ',trim(file_path),' does not exist' - endif + endif end function @@ -159,47 +217,128 @@ logical function symbol_search(line,symbol) integer i character(len) line ; character(1) symbol ; symbol_search = .false. do i=1,len - if (line(i:i)==symbol) symbol_search = .true. - end do + if (line(i:i)==symbol) then + symbol_search = .true. + exit + endif + end do end function logical function word_search(line,goal) !=- of only the first word - ? character(len) word,line ; character(*) goal ; word_search = .false. read(line,*) word - if (word==goal) word_search = .true. + if (word==goal) word_search = .true. + end function + + logical function word_search_array(array,goal) !=- of only the first word - ? + character(len) word,array(:) ; character(*) goal ; integer i + word_search_array = .false. + do i=1,size(array) + read(array(i),*) word + if (word==goal) word_search_array = .true. + enddo end function subroutine shell_MD_Tree + integer i do i=1,size(folders(:)) - select case (operation_system) - case(1) !=- Windows - folders(i) = backslashfix(folders(i)) - if (folders(i).ne.'') call system( 'MD ' // trim(folders(i)) ) - case(2) !=- Linux - folders(i) = slashfix(folders(i)) - if (folders(i).ne.'') call system( 'mkdir -p ' // trim(folders(i)) ) - end select + if (folders(i).ne.'') call create_folder(folders(i)) enddo end subroutine shell_MD_Tree + subroutine move_file(whence,here) + character(len) whence , here + select case (operating_system) + case(1) !=- Windows + whence = backslashfix(whence) + here = backslashfix(here) + if ( whence.ne.'' .and. here.ne.'' .and. whence.ne.'NaN' .and. here.ne.'NaN' ) & + call system( 'move ' // trim(whence) // ' ' // trim(here) // ' >> log.log ' ) + case(2) !=- Linux + whence = slashfix(whence) + here = slashfix(here) + if ( whence.ne.'' .and. here.ne.'' .and. whence.ne.'NaN' .and. here.ne.'NaN' ) & + call system( 'mv ' // trim(whence) // ' ' // trim(here) // ' >> log.log ' ) + end select + end subroutine + + + subroutine copy_file(whence,here) + character(len) whence , here + select case (operating_system) + case(1) !=- Windows + whence = backslashfix(whence) + here = backslashfix(here) + if ( whence.ne.'' .and. here.ne.'' .and. whence.ne.'NaN' .and. here.ne.'NaN' ) & + call system( 'copy ' // trim(whence) // ' ' // trim(here) // ' >> log.log ' ) + case(2) !=- Linux + whence = slashfix(whence) + here = slashfix(here) + if ( whence.ne.'' .and. here.ne.'' .and. whence.ne.'NaN' .and. here.ne.'NaN' ) & + call system( 'cp ' // trim(whence) // ' ' // trim(here) // ' >> log.log ' ) + end select + end subroutine + + + subroutine create_folder(path_folder) + character(len) path_folder + select case (operating_system) + case(1) !=- Windows + path_folder = backslashfix(path_folder) + if (path_folder.ne.'NaN' .and. path_folder.ne.'') & + call system( 'If Not Exist ' // trim(path_folder) // & + ' MD ' // trim(path_folder) // ' >> log.log ' ) + case(2) !=- Linux + path_folder = slashfix(path_folder) + if (path_folder.ne.'NaN' .and. path_folder.ne.'') & + call system( 'mkdir -p ' // trim(path_folder) // ' >> log.log ' ) + end select + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + character(len) function slashfix(path) integer i character(len) path; slashfix = ''; slashfix = path do i=len,1,-1 ; if (path(i:i)=='\') slashfix(i:i)='/' ; enddo end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + character(len) function backslashfix(path) integer i character(len) path; backslashfix = ''; backslashfix = path do i=len,1,-1 ; if (path(i:i)=='/') backslashfix(i:i)='\' ; enddo end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + integer function count_noNaN(array) + real(8) array(:) + count_noNaN = 0 + do i=1,size(array) + if ( array(i) .ne. NaN ) count_noNaN = 1 + count_noNaN + enddo + if ( count_noNaN==0 .and. all_error_reporst ) write(*,*) 'count_noNaN: 0' + end function count_noNaN + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + integer function count_noNaN_text(array) + character(len) array(:) + count_noNaN_text = 0 + do i=1,size(array) + if ( array(i) .ne. 'NaN' ) count_noNaN_text = 1 + count_noNaN_text + enddo + if ( count_noNaN_text==0 .and. all_error_reporst ) write(*,*) 'count_noNaN_text: 0' + end function count_noNaN_text + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! - end module !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/main.f95 b/main.f95 index b75f16a..3b95834 100644 --- a/main.f95 +++ b/main.f95 @@ -1,7 +1,6 @@ !=- fortran-libraries !=- © Stanislav Shirokov, 2014-2020 -!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 program COSMOD use COSMOD_citadel @@ -9,6 +8,5 @@ program COSMOD call commands -!=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 end !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 diff --git a/math.f95 b/math.f95 index eb58992..bff56af 100644 --- a/math.f95 +++ b/math.f95 @@ -4,44 +4,469 @@ !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 module math use global + use GNUplot + use Cosmology - integer,parameter :: N_medians = 1d2 + logical :: approx_save_datafiles = .false. - integer :: N_grid_LT = 2d2, N_grid_LT_a = 2d2, N_grid_LT_b = 2d1, N_points + character(len) :: columns_for_approx = 'NaN' , & !=- x y | x y dy | x y dx dy + polylog_a_data = 'polylog_a.dat' , & + approx_datafile = 'data_approx.dat' , & + approx_coefficients = 'approx_coefficients.dat' - real(8) :: RightLimit = 1d150, LeftLimit = -1d150, pi = 4*datan(1d0), & + integer,parameter :: N_medians = 1d2 , approx_max_order = 10 , N_temp_approx = 1d3 + + real(8),parameter :: pi = 4*datan(1d0) + + integer :: N_grid_LT = 2d2, N_grid_LT_a = 2d2, N_grid_LT_b = 2d1, N_points , & + data_column_x , data_column_y , data_column_dx , data_column_dy, & + + approximating_step , step_N_temp , approx_order = 2 , & + approx_minimum_order = 99 + + real(8) :: RightLimit = 1d150, LeftLimit = -1d150, & a_error_lower, b_error_lower, a_error_upper, b_error_upper, a_start = 1d-4, b_start = 1d-4, a_max = 1d1, b_max = 1d0, & - XYdXdY(4,N_GRB), mx, my, a_step, b_step, X2=1d150, va, vb, weights, Sx, Sy, & + mx, my, a_step, b_step, X2=1d150, va, vb, weights, Sx, Sy, & va_error_lower,vb_error_lower , va_error_upper, vb_error_upper, medar(4,N_medians), & add_median_left_border, add_median_right_border, x_step,x_step_start,x_step_final, & add_median_x,add_median_y,add_median_dx,add_median_dy, log_medians_parameter = 0, & - add_median2_x,add_median2_y,add_median2_dx,add_median2_dy + add_median2_x,add_median2_y,add_median2_dx,add_median2_dy, & + + polylog_a(approx_max_order+1) = 0d0 , polylog_a_borders(2,approx_max_order+1) = 1d0 , approx_scaling_x = 1d2 , & + extremums_polylog_a( approx_max_order+1 , approx_max_order+1 ) , & + extremums_chi2(approx_max_order+1) = NaN, linear_approx(2) = 0d0 , minimum_chi2 = NaN, & + temp_a(approx_max_order+1) = 0d0 , approx_minimum_chi2 = NaN , approx_minimum_polylog_a(approx_max_order+1) = NaN + + real(8),allocatable,dimension(:,:) :: XYdXdY !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1 contains + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine math_default + + columns_for_approx = 'NaN' + polylog_a_data = 'polylog_a.dat' + approx_datafile = 'data_approx.dat' + approx_coefficients = 'approx_coefficients.dat' + approx_minimum_order = 99 + approx_minimum_chi2 = NaN + approx_minimum_polylog_a(:) = NaN + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine approx_default + + polylog_a_data = 'polylog_a.dat' + approx_datafile = 'data_approx.dat' + approx_coefficients = 'approx_coefficients.dat' + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine approx_columns_checking + integer :: x = 0 , y = 0 , dx = 0 , dy = 0 + + if (columns_for_approx=='NaN') then + x=1 ; y=2 + else + + read(columns_for_approx,*,end=11,err=11) x , y + read(columns_for_approx,*,end=11) x , y , dy + read(columns_for_approx,*,end=11) x , y , dy , dx + 11 continue ; if ( x==0 .or. y==0 ) write(*,*) ' scripting_error: approx_columns_checking: x=0 or y=0' + + endif + data_column_x = x + data_column_y = y + data_column_dx = dx + data_column_dy = dy + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine read_approx_catalog(data_file) + character(len) :: data_file , line , catalog(200) = '0' + integer i + + unit_1 = random_unit() + open(unit_1,file=data_file,status='old') + + i=1 + read(unit_1,'(A4096)',end=22) line !=- 1-st line reading (head of data table) + do + read(unit_1,'(A4096)',end=22) line + read(line,*,end=11,err=11) catalog + 11 continue + + read(catalog( data_column_x ),*) XYdXdY( 1 , i ) + read(catalog( data_column_y ),*) XYdXdY( 2 , i ) + if (data_column_dx/=0) read(catalog( data_column_dx ),*) XYdXdY( 3 , i ) + if (data_column_dy/=0) read(catalog( data_column_dy ),*) XYdXdY( 4 , i ) + i=1+i + + enddo ; 22 continue + close(unit_1) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + real(8) function polylog( a , x , order ) + integer i,order + real(8) x,a(order+1) + + polylog = 0d0 + do i=0,order + polylog = polylog + a(1+i) * log10 (x) ** i + !write(*,*) temp_a(1:2) , temp_a(2) * log10 (x) + enddo + + polylog = polylog + 5d0*log10 (x) !+ 25d0 !- 1d1 + + end function + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + real(8) function chi2_polylog(order) + integer i,order,N + real(8) polylog_model + + N = size(XYdXdY(1,:)) + + ! XYdXdY( 4 , : ) = 1d0 !=- test + + chi2_polylog = 0d0 + do i=1,N + if (XYdXdY( 4 , i )==0) XYdXdY( 4 , i ) =1d0 + polylog_model = polylog( temp_a , XYdXdY( 1 , i ) , order ) + chi2_polylog = chi2_polylog + ( XYdXdY( 2 , i ) - polylog_model ) ** 2 / & + XYdXdY( 4 , i )**2 / dabs(polylog_model) / N + enddo + + end function + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine write_polylog_a( order ) + integer order ; real(8) chi2 + + chi2 = chi2_polylog( order ) + chi2 = log10(chi2) + + unit_2 = random_unit() + open(unit_2,file=polylog_a_data,status='old',position='append') + + theformat = '(' // trim(adjustl(inttostr(approx_max_order))) // '(E20.8))' + write(unit_2,theformat) chi2 , temp_a(1:order+1) + + if ( chi2 < minimum_chi2 ) then + minimum_chi2 = chi2 + polylog_a(:) = temp_a(:) + endif + + close(unit_2) + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine create_data_polylog_a_temp + unit_2 = random_unit() + open(unit_2,file=polylog_a_data,status='replace') + close(unit_2) + end subroutine + + subroutine write_data_approx( order ) + integer i,order + real(8) :: temp_a(7) = (/ 33.158168604381189 , & + 6.2594070261452195E-002 , & + 2.7994852249554239E-002 , & + 2.7962159075643332E-002 , & + 1.8916688159884384E-002 , & + 3.0028150691020179E-002 , & + -1.2242181029368278E-002 /) + + !=- polylog_a(1:4) = (/ 44.1 , 6.28 , 0.57 , 0.08 /) + + if (approx_save_datafiles) then + approx_datafile = trim(file_path(approx_datafile)) // trim(adjustl( inttostr(order) )) // '_' // & + trim(name_ext(approx_datafile)) + approx_coefficients = trim(file_path(approx_coefficients)) // trim(adjustl( inttostr(order) )) // '_' // & + trim(name_ext(approx_coefficients)) + end if + + unit_3 = random_unit() + open(unit_3,file=approx_datafile,status='replace') + do i=1,size(XYdXdY(1,:)) + write(unit_3,*) XYdXdY(1,i)/approx_scaling_x, XYdXdY(2:4,i) , polylog( polylog_a , XYdXdY( 1 , i ) , order ) , & + linear_approx(1) + log10(XYdXdY(1,i))*(linear_approx(2)) , dm_LCDM( XYdXdY(1,i)/approx_scaling_x , 7d1 ), & + polylog( temp_a , XYdXdY( 1 , i ) , size(temp_a)-1 ) , 43d0 + 5d0*log10( XYdXdY(1,i)/approx_scaling_x ) + enddo + close(unit_3) + + unit_4 = random_unit() + open(unit_4,file=approx_coefficients,status='replace') + write(unit_4,*) order , 1d1**minimum_chi2 , polylog_a(1:order+1) + close(unit_4) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine approx_polylog(data_file) + character(len) data_file + real(8) :: shift(approx_max_order+1) = NaN , center , target_shift , accuracy = 1d-4 + integer i,N,order , condition + + N = file_volume(data_file) - 1 + allocate ( XYdXdY( 4 , N ) ) ; XYdXdY(:,:)=1d0 + + call approx_columns_checking + call read_approx_catalog(data_file) + + XYdXdY(1,:) = approx_scaling_x * XYdXdY(1,:) + + call logtrend_x(XYdXdY(1:2,:),linear_approx(2),linear_approx(1),0d0*approx_scaling_x,0.1*approx_scaling_x) + !=- write(*,*) ' LA: ' , linear_approx(1) + 25 , linear_approx(2) - 5 + + if (approx_order<0) approx_order=0 + if (approx_order>approx_max_order) approx_order=approx_max_order + do j=1,approx_order + + call approx_default + + order = j + minimum_chi2 = NaN + approximating_step = 1 + condition = 0 + + polylog_a_borders(1,1) = 0.9 * ( linear_approx(1) ) !=- + 25 + polylog_a_borders(2,1) = 1.1 * ( linear_approx(1) ) !=- + 25 + polylog_a_borders(1,2:approx_max_order) = -5d0 + polylog_a_borders(2,2:approx_max_order) = 5d0 + + do while ( condition < 3 ) + + call create_data_polylog_a_temp + + do i = 1,N_temp_approx + step_N_temp = i + call random_polylog_a(N_temp_approx) + call write_polylog_a( order ) + enddo + + do i=1,order+1 + center = polylog_a_borders(1,i) + (polylog_a_borders(2,i) - polylog_a_borders(1,i))*0.5d0 + target_shift = ( polylog_a_borders(2,i) - center ) * 0.9d0 + shift(i) = polylog_a(i) - center + + if ( shift(i) < target_shift ) then + polylog_a_borders(1,i) = polylog_a(i) - target_shift + polylog_a_borders(2,i) = polylog_a(i) + target_shift + else + polylog_a_borders(:,i) = polylog_a_borders(:,i) + shift(i) + end if + enddo + + call plot_polylog_a( order ) + + extremums_polylog_a(j,:) = polylog_a(:) + extremums_chi2(j) = minimum_chi2 + approximating_step = 1 + approximating_step + + if ( minimum_chi2 < approx_minimum_chi2 ) then + approx_minimum_chi2 = minimum_chi2 + approx_minimum_polylog_a(:) = polylog_a(:) + approx_minimum_order = j + endif + + if (sum(dabs(shift(1:order+1))) < accuracy*(order+1)) condition = 1 + condition + + enddo + call write_data_approx( order ) + call plot_data_approx( data_file , order ) + end do + write(*,*) ' math_approx: the best order:', approx_minimum_order + write(*,*) ' math_approx: the best coefficients:', & + approx_minimum_polylog_a(1) , approx_minimum_polylog_a(2:approx_minimum_order+1) + write(*,*) ' math_approx: the best chi2:', 1d1**approx_minimum_chi2 , approx_minimum_chi2 + deallocate(XYdXdY) + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine plot_polylog_a( order ) + integer order + call clear_plot_fields + + graph_name = trim(adjustl(inttostr(order))) // '_' // trim(adjustl(inttostr(approximating_step))) + + pngterm = 'set term pngcairo enhanced font "Verdana,20" size 1024, 1024' + epsterm = 'set term postscript portrait enhanced color font "Verdana,10" size 6, 6' + + GNUfields (legend) = 'off' + + GNUfields (add1) = 'set palette rgb 7,5,15' + GNUfields (add2) = 'set view map' + GNUfields (add3) = 'set cblabel "log {/Symbol C}^2"' + GNUfields (add4) = 'set format cb "%2.0t{/Symbol \327}10^{%L}"' + + GNUfields (title) = '{/Symbol C}^2-phase plane a_0--a_1, approximating step ' // & + trim(adjustl(inttostr(approximating_step))) // ' (CosMod v1.0)' + GNUfields (xlabel) = 'a_0' + GNUfields (ylabel) = 'a_1' + + GNUfields (plot1) = 'splot "' // trim((polylog_a_data)) // & + '" u 2:3:1 w p lw 2 pt 7 ps 0.5 dt 1 palette' + + GNUfields(extention_out_figure)='eps' ; call plot(polylog_a_data) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine plot_data_approx( source_data , order ) + integer order + character(len) source_data + character(2) :: delta_column = '7' + + call clear_plot_fields + + graph_name = trim(adjustl(inttostr(order))) // '_' // trim(name(source_data)) + + pngterm = 'set term pngcairo enhanced font "Verdana,20" size 1400, 1050' + epsterm = 'set term postscript enhanced color font "Verdana,14" size 8.5, 6.3' + + GNUfields (legend) = 'right bottom' + + GNUfields (title) = ' approximation of the ' // trim(name(source_data)) // & + ' data (CosMod v1.0)' + GNUfields (xlabel) = 'z' + GNUfields (ylabel) = '{/Symbol m}' + GNUfields (logscale) = 'set logscale x' + + GNUfields (ls1) = 'set linestyle 1 lw 1 pt 7 ps 0.7 lt rgb "blue"' + GNUfields (ls2) = 'set linestyle 2 lw 2 pt 7 ps 0.7 lt rgb "black"' + GNUfields (ls3) = 'set linestyle 3 lw 4 pt 7 ps 0.7 lt rgb "yellow"' + GNUfields (ls4) = 'set linestyle 4 lw 3 pt 7 ps 0.7 lt rgb "red"' + GNUfields (ls5) = 'set linestyle 5 lw 3 pt 7 ps 0.7 lt rgb "dark-cyan"' + + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:2:4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:5 w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:6 w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:7 w l ls 4' + + GNUfields (title1) = 'the ' // trim(name(source_data)) // ' data' + GNUfields (title2) = 'polylog approximation of power ' // trim( adjustl( inttostr(order) ) ) + GNUfields (title3) = 'first linear approximation (FLA)' + GNUfields (title4) = 'LCDM (70,0.7)' + + + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) + + graph_name = 'delta_LCDM_' // trim(adjustl(inttostr(order))) // '_' // trim(name(source_data)) + + GNUfields (ylabel) = '{/Symbol D}{/Symbol m} = {/Symbol m} - {/Symbol m}_{LCDM(70,0.7)}' + + delta_column = '7' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:($2-$' // delta_column // '):4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:($5-$' // delta_column // ') w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:($6-$' // delta_column // ') w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:($7-$' // delta_column // ') w l ls 4' + GNUfields (plot5) = '"' // trim((approx_datafile)) // '" u 1:($8-$' // delta_column // ') w l ls 5' + + GNUfields (title5) = 'testing model' + + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) + + graph_name = '5log_x_' // trim(adjustl(inttostr(order))) // '_' // trim(name(source_data)) + + GNUfields (ylabel) = '{/Symbol D}{/Symbol m} = {/Symbol m} - {/Symbol m}_{55+5log x}' + + delta_column = '8' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:($2-$' // delta_column // '):4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:($5-$' // delta_column // ') w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:($6-$' // delta_column // ') w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:($7-$' // delta_column // ') w l ls 4' + GNUfields (plot5) = '"' // trim((approx_datafile)) // '" u 1:($8-$' // delta_column // ') w l ls 5' + + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) + + graph_name = 'delta_linear_' // trim(adjustl(inttostr(order))) // '_' // trim(name(source_data)) + + GNUfields (ylabel) = '{/Symbol D}{/Symbol m} = {/Symbol m} - {/Symbol m}_{FLA}' + GNUfields (legend) = 'left top' + + delta_column = '6' + GNUfields (plot1) = '"' // trim((approx_datafile)) // '" u 1:($2-$' // delta_column // '):4 w yerrorb ls 1' + GNUfields (plot2) = '"' // trim((approx_datafile)) // '" u 1:($5-$' // delta_column // ') w l ls 3' + GNUfields (plot3) = '"' // trim((approx_datafile)) // '" u 1:($6-$' // delta_column // ') w l ls 2' + GNUfields (plot4) = '"' // trim((approx_datafile)) // '" u 1:($7-$' // delta_column // ') w l ls 4' + GNUfields (plot5) = '"' // trim((approx_datafile)) // '" u 1:($8-$' // delta_column // ') w l ls 5' + + GNUfields(extention_out_figure)='png' ; call plot(approx_datafile) + + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine random_polylog_a(N_temp) + real(8) x , stretching , center + integer i , N_temp + + do i=1,approx_max_order + + call random_number(x) + x=x-0.5d0 + + center = polylog_a_borders(1,i) + (polylog_a_borders(2,i) - polylog_a_borders(1,i))*0.5d0 + stretching = x*step_N_temp**2 * ( polylog_a_borders(2,i) - polylog_a_borders(1,i) ) / N_temp**2 + + + temp_a(i) = center + stretching + !temp_a(i) = polylog_a_borders(1,i) + a * ( polylog_a_borders(2,i) - polylog_a_borders(1,i) ) + + end do + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real function distance(x,y,z) real(8) x,y,z distance = ( x**2 + y**2 + z**2 )**0.5d0 end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + integer function sign_real(x) ; real(8) x; sign_real=1 ; if (x<0) sign_real=-1 end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function mydexp(x) ; real(8) x; n=1 ; mydexp = ( 1d0 + x/1d5 )**1d5 end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function factorial(n) ; integer n; factorial=1d0; do while (n>0); factorial=factorial*n; n=n-1; enddo end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + real(8) function mydsin(x) ; real(8) x; n=1; mydsin=x; do; mydsin=mydsin+(-1d0)**n * x**(2*n+1) / factorial (2*n+1) n=n+1 ; if (x**(2*n+1) / factorial (2*n+1) < 1d-19) exit ; end do end function + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine logtrend(XY,a,b,Lb,Rb) real(8) XY(2,N_GRB),G(2,N_GRB),a,b,sumx,sumx2,sumy,sumxy,Lb,Rb; G(:,:)=0d0; forall (i=1:N_GRB,j=1:2, XY(j,i)>0) G(j,i)=log10(XY(j,i)) @@ -51,6 +476,27 @@ subroutine logtrend(XY,a,b,Lb,Rb) if (k.ne.0) then ; a=(k*sumxy-sumx*sumy)/(k*sumx2-sumx**2) ; b=(sumy-a*sumx)/k ; endif end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + + subroutine logtrend_x(XY,a,b,Lb,Rb) + integer N,i,j,k + real(8) XY(:,:) ,a,b,sumx,sumx2,sumy,sumxy,Lb,Rb + real(8),allocatable,dimension(:,:) :: G + + N = size(XY(2,:)) + allocate(G(2,N)) + G(:,:)=XY(:,:); forall (i=1:N, XY(1,i)>0) G(1,i)=log10(XY(1,i)) + + k=0;sumx=0d0;sumx2=0d0;sumy=0d0;sumxy=0d0 + do j=1,N; if (G(2,j).ne.NaN .and. G(1,j).le.log10(Rb) .and. G(1,j).ge.log10(Lb)) then! write(*,*) N,G(1,j),G(2,j),XY(1,j),XY(2,j) + sumx=sumx+G(1,j) ; sumx2=sumx2+G(1,j)**2d0 ; sumy=sumy+G(2,j) ; sumxy=sumxy+G(1,j)*G(2,j) ; k=k+1 ; endif ; enddo + if (k.ne.0) then ; a=(k*sumxy-sumx*sumy)/(k*sumx2-sumx**2) ; b=(sumy-a*sumx)/k ; endif + + deallocate(G) + end subroutine + + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine logtrendX2(Array) real(8) Array(:,:) !=- XY(1,i) - x, XY(3,i) - dx @@ -66,7 +512,7 @@ subroutine logtrendX2(Array) mx = sum(XYdXdY(1,:))/N_points my = sum(XYdXdY(2,:))/N_points -if ( abs(mx) < 1 ) write(*,*) 'bad solution' + if ( abs(mx) < 1 ) write(*,*) 'bad solution' forall (i=1:N, Array(1,i)>0) XYdXdY(1,i) = XYdXdY(1,i) - mx forall (i=1:N, Array(2,i)>0) XYdXdY(2,i) = XYdXdY(2,i) - my @@ -121,6 +567,8 @@ subroutine logtrendX2(Array) end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine SORT2(Array) real(8) Array(:,:),buf(2) do i=1,size(Array(1,:)) @@ -131,6 +579,8 @@ subroutine SORT2(Array) enddo end subroutine + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine SORT(Array) real(8) Array(:),abuf do i=1,size(Array(:)) @@ -141,7 +591,7 @@ subroutine SORT(Array) enddo end subroutine - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine compute_medians(input_array,x_start_in,x_final_in,dx) real(8), intent(in) :: input_array(:,:),dx,x_start_in,x_final_in @@ -253,7 +703,7 @@ subroutine compute_medians(input_array,x_start_in,x_final_in,dx) deallocate(XY) end subroutine - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! real(8) function median(Array) real(8) Array(:) @@ -276,7 +726,7 @@ real(8) function median(Array) deallocate(Barray) end function - + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! subroutine EclCoord(RA,Dec,l2,b2) != from galactic to ecliptical ==================================1 implicit none ;real(8) RA,RA2,Dec,l,l2,b,b2,RAl,Decl,l0,pi,r,d @@ -287,6 +737,8 @@ subroutine EclCoord(RA,Dec,l2,b2) != from galactic to ecliptical =============== Dec=Dec*180d0/pi;if (RA2.le.0) RA=pi-RA;RA=RA*180d0/pi-RAl + 25.71667480468742d0;if (RA.le.0) RA=RA+360 end subroutine EclCoord !==========================================================================1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine GalCoord(RA,Dec,l,b) != from ecliptical to galactic ====================================1 implicit none ;real(8) RA,Dec,l,l2,b,RAl,Decl,l0,pi,r,d pi=4*datan(1d0);l0=32.93192+90;RAl=192.858333/180d0*pi;Decl=27.128333/180d0*pi @@ -295,6 +747,8 @@ subroutine GalCoord(RA,Dec,l,b) != from ecliptical to galactic ================= if (l2.le.0) l=pi-l;l=l0-l*180d0/pi;if (l.le.0) l=l+360;RA=r;Dec=d;b=b*180d0/pi end subroutine GalCoord !==========================================================================1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine SphCoord(x,y,z,b,l,r) != from x,y,z to a,b,r ===========================================1 real(8) x,y,z,l,b,r !=- b -> l , l -> b, =1 r=dsqrt(x**2+y**2+z**2);if (r==0) then ; l=0; b=0; goto 11 ; else ; l=datan(z/dsqrt(x**2+y**2)) @@ -303,6 +757,8 @@ subroutine SphCoord(x,y,z,b,l,r) != from x,y,z to a,b,r ======================== l=l*180d0/pi ; b=b*180d0/pi end subroutine !===================================================================================1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine SphCoord_radians(x,y,z,l,b,r) != from x,y,z to a,b,r ===========================================1 real(8) x,y,z,l,b,r r=dsqrt(x**2+y**2+z**2);if (r==0) then ; l=0; b=0; goto 11 ; else ; l=datan(z/dsqrt(x**2+y**2)) @@ -310,9 +766,13 @@ subroutine SphCoord_radians(x,y,z,l,b,r) != from x,y,z to a,b,r ================ if (x==0.and.y.gt.0) b=pi/2;if (x==0.and.y.lt.0) b=3d0*pi/2 ; endif; 11 continue end subroutine !===================================================================================1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + subroutine DescCoord(a,b,r,x,y,z) != from a,b,r to x,y,z ==========================================1 implicit none ; real(8) x,y,z,a,b,r ; x=r*dcos(a)*dcos(b) ; y=r*dcos(a)*dsin(b);z=r*dsin(a) end subroutine !===================================================================================1 + !=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=! + end module !=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- truncated=136-=1