In [1]:

    

cd /astrolab21/gabriel/mesa/r8845/make_co_wd_r8845/


    

    




/astrolab21/gabriel/mesa/r8845/make_co_wd_r8845

In [3]:

    

%matplotlib inline
import mesa as ms
import numpy as np
import matplotlib
matplotlib.rcParams['figure.figsize'] = (14.0, 7.0) #(14.0, 10.0)
import matplotlib.pyplot as plt


    

In [4]:

    

# Lendo os dados
folder = 'LOGS/'
hist = ms.history_data(folder)
teff = hist.get('log_Teff')
lum = hist.get('log_L')
models = hist.get('model_number')

# Plot do HRD
plt.plot(teff, lum)
plt.gca().invert_xaxis()
plt.ylabel(r'$\log L$')
plt.xlabel(r'$\log Teff$')


    

    




No history.datasa file found, create new one from history.data
 reading ...100% 







    
Out[4]:





<matplotlib.text.Text at 0x7fb88f3dc890>






    

In [7]:

    

# Zoom no HRD entre log Teff 3.5, 3.7
plt.plot(teff, lum)
plt.gca().invert_xaxis()
plt.ylabel(r'$\log L$')
plt.xlabel(r'$\log Teff$')
plt.xlim(3.7, 3.5)
plt.ylim(3, 4)


    

    
Out[7]:





(3, 4)






    

In [38]:

    

# Perfil quimico no ultimo modelo
prof = ms.mesa_profile(folder, models[-1], num_type='nearest_model')

q = prof.get('q')
logq = - prof.get('logxq')
mass = hist.get('star_mass')
tfinal = 10**teff[-1]

h1 = prof.get('h1')
#he3 = prof.get('he3')
he4 = prof.get('he4')
c12 = prof.get('c12')
#c13 = prof.get('c13')
#n13 = prof.get('n13')
n14 = prof.get('n14')
#n15 = prof.get('n15')
o16 = prof.get('o16')
#o18 = prof.get('o18')
ne20 = prof.get('ne20')
#ne22 = prof.get('ne22')
mg24 = prof.get('mg24')

plt.plot(logq, he4, color='darkolivegreen', lw=1.5, label=r'$He_4$')
plt.plot(logq, c12, color='k', lw=1.5, label=r'$C_{12}$')
plt.plot(logq, n14, color='r', lw=1.5, label=r'$N_{14}$')
plt.plot(logq, o16, color='g', lw=1.5, label=r'$O_{16}$')
plt.plot(logq, ne20, color='y', lw=1.5, label=r'$Ne_{20}$')
plt.plot(logq, mg24, color='m', lw=1.5, label=r'$Mg_{24}$')
plt.plot(logq, h1, color='b', lw=1.5, label=r'$H_1$')
plt.xlabel(r'$- \log (1 - q)$')
plt.ylabel('mass fraction')
plt.title('Perfil Quimico ' + ' - model n ' + str(models[-1]+1))
plt.suptitle(r'$M_f = $ ' + str(mass[-1]) + r'-  $T_{eff} = $ ' + str(tfinal))
plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
plt.xlim(0, 20)
plt.savefig('make_co_final_abun.png')


    

    




66 in profiles.index file ...
Found and load nearest profile for cycle 3437
reading LOGS//profile66.data ...
 reading ...100% 

Closing profile tool ...






    

In [26]:

    

# star mass vs model
plt.plot(models, mass)
print 'last model = ', models[-1]


    

    




last model =  3437.0
(232,)






    

In [34]:

    

# star_teff vs model
teffK = 10**teff
plt.plot(models, teffK)


    

    
Out[34]:





[<matplotlib.lines.Line2D at 0x7fb886f20390>]






    

In [35]:

    

cat src/run_star_extras.f


    

    




! ***********************************************************************
!
!   Copyright (C) 2010  Bill Paxton
!
!   this file is part of mesa.
!
!   mesa is free software; you can redistribute it and/or modify
!   it under the terms of the gnu general library public license as published
!   by the free software foundation; either version 2 of the license, or
!   (at your option) any later version.
!
!   mesa is distributed in the hope that it will be useful, 
!   but without any warranty; without even the implied warranty of
!   merchantability or fitness for a particular purpose.  see the
!   gnu library general public license for more details.
!
!   you should have received a copy of the gnu library general public license
!   along with this software; if not, write to the free software
!   foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************
 
      module run_star_extras

      use star_lib
      use star_def
      use const_def
      
      implicit none
      
      integer :: time0, time1, clock_rate
      double precision, parameter :: expected_runtime = 6.8 ! minutes

      
      contains
      
      
      subroutine extras_controls(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         s% extras_startup => extras_startup
         s% extras_check_model => extras_check_model
         s% extras_finish_step => extras_finish_step
         s% extras_after_evolve => extras_after_evolve
         s% how_many_extra_history_columns => how_many_extra_history_columns
         s% data_for_extra_history_columns => data_for_extra_history_columns
         s% how_many_extra_profile_columns => how_many_extra_profile_columns
         s% data_for_extra_profile_columns => data_for_extra_profile_columns  
      end subroutine extras_controls
      
      
      integer function extras_startup(id, restart, ierr)
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_startup = 0
         call system_clock(time0,clock_rate)
         if (.not. restart) then
            call alloc_extra_info(s)
         else ! it is a restart
            call unpack_extra_info(s)
         end if
      end function extras_startup
      
      
      subroutine extras_after_evolve(id, id_extra, ierr)
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         double precision :: dt
         character (len=strlen) :: test
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         call system_clock(time1,clock_rate)
         dt = dble(time1 - time0) / clock_rate / 60
         call GET_ENVIRONMENT_VARIABLE( &
            "MESA_TEST_SUITE_CHECK_RUNTIME", test, status=ierr, trim_name=.true.)
         if (ierr == 0 .and. trim(test) == 'true' .and. dt > 1.5*expected_runtime) then
            write(*,'(/,a70,2f12.1,99i10/)') &
               'failed: EXCESSIVE runtime, prev time, retries, backups, steps', &
               dt, expected_runtime, s% num_retries, s% num_backups, s% model_number
         else
            write(*,'(/,a50,2f12.1,99i10/)') 'runtime, prev time, retries, backups, steps', &
               dt, expected_runtime, s% num_retries, s% num_backups, s% model_number
         end if
         ierr = 0
      end subroutine extras_after_evolve
      

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         
         real(dp), parameter :: Blocker_scaling_factor_after_TP = 5d0
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         include 'formats'
         
         extras_check_model = keep_going
         
         if (abs(s% Blocker_scaling_factor - Blocker_scaling_factor_after_TP) < 1d-8) return
         
         if (s% center_he4 < 1d-4 .and. &
               any(s% burn_he_conv_region(1:s% num_conv_boundaries)) .and. &
               s% he_core_mass - s% c_core_mass <= s% TP_he_shell_max) then
            !write(*,1) 'set Blocker_scaling_factor = Blocker_scaling_factor_after_TP', Blocker_scaling_factor_after_TP
            !s% Blocker_scaling_factor = Blocker_scaling_factor_after_TP
            extras_check_model = terminate
         end if
         
      end function extras_check_model


      integer function how_many_extra_history_columns(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 0
      end function how_many_extra_history_columns
      
      
      subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine data_for_extra_history_columns

      
      integer function how_many_extra_profile_columns(id, id_extra)
         use star_def, only: star_info
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_columns = 0
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(id, id_extra, n, nz, names, vals, ierr)
         use star_def, only: star_info, maxlen_profile_column_name
         use const_def, only: dp
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine data_for_extra_profile_columns
      

      ! returns either keep_going or terminate.
      integer function extras_finish_step(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_finish_step = keep_going
         call store_extra_info(s)
      end function extras_finish_step
      
      
      ! routines for saving and restoring extra data so can do restarts
         
         ! put these defs at the top and delete from the following routines
         !integer, parameter :: extra_info_alloc = 1
         !integer, parameter :: extra_info_get = 2
         !integer, parameter :: extra_info_put = 3
      
      
      subroutine alloc_extra_info(s)
         integer, parameter :: extra_info_alloc = 1
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_alloc)
      end subroutine alloc_extra_info
      
      
      subroutine unpack_extra_info(s)
         integer, parameter :: extra_info_get = 2
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_get)
      end subroutine unpack_extra_info
      
      
      subroutine store_extra_info(s)
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_put)
      end subroutine store_extra_info
      
      
      subroutine move_extra_info(s,op)
         integer, parameter :: extra_info_alloc = 1
         integer, parameter :: extra_info_get = 2
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         integer, intent(in) :: op
         
         integer :: i, j, num_ints, num_dbls, ierr
         
         i = 0
         ! call move_int or move_flg    
         num_ints = i
         
         i = 0
         ! call move_dbl       
         
         num_dbls = i
         
         if (op /= extra_info_alloc) return
         if (num_ints == 0 .and. num_dbls == 0) return
         
         ierr = 0
         call star_alloc_extras(s% id, num_ints, num_dbls, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in star_alloc_extras'
            write(*,*) 'alloc_extras num_ints', num_ints
            write(*,*) 'alloc_extras num_dbls', num_dbls
            stop 1
         end if
         
         contains
         
         subroutine move_dbl(dbl)
            double precision :: dbl
            i = i+1
            select case (op)
            case (extra_info_get)
               dbl = s% extra_work(i)
            case (extra_info_put)
               s% extra_work(i) = dbl
            end select
         end subroutine move_dbl
         
         subroutine move_int(int)
            integer :: int
            i = i+1
            select case (op)
            case (extra_info_get)
               int = s% extra_iwork(i)
            case (extra_info_put)
               s% extra_iwork(i) = int
            end select
         end subroutine move_int
         
         subroutine move_flg(flg)
            logical :: flg
            i = i+1
            select case (op)
            case (extra_info_get)
               flg = (s% extra_iwork(i) /= 0)
            case (extra_info_put)
               if (flg) then
                  s% extra_iwork(i) = 1
               else
                  s% extra_iwork(i) = 0
               end if
            end select
         end subroutine move_flg
      
      end subroutine move_extra_info

      end module run_star_extras
      

In [36]:

    

ls LOGS/


    

    




history.data    profile21.data  profile34.data  profile47.data  profile5.data
history.datasa  profile22.data  profile35.data  profile48.data  profile60.data
profile10.data  profile23.data  profile36.data  profile49.data  profile61.data
profile11.data  profile24.data  profile37.data  profile4.data   profile62.data
profile12.data  profile25.data  profile38.data  profile50.data  profile63.data
profile13.data  profile26.data  profile39.data  profile51.data  profile64.data
profile14.data  profile27.data  profile3.data   profile52.data  profile65.data
profile15.data  profile28.data  profile40.data  profile53.data  profile66.data
profile16.data  profile29.data  profile41.data  profile54.data  profile6.data
profile17.data  profile2.data   profile42.data  profile55.data  profile7.data
profile18.data  profile30.data  profile43.data  profile56.data  profile8.data
profile19.data  profile31.data  profile44.data  profile57.data  profile9.data
profile1.data   profile32.data  profile45.data  profile58.data  profiles.index
profile20.data  profile33.data  profile46.data  profile59.data

In [ ]: