anchor.diff:       &        halo_stars_fraction,         burst_mrep,
        &        normal_mrep,                 thick_disk_max_sfr,
        &        thin_disk_birth_init_time,   thick_disk_age,
  -     &        thick_disk_max_sfr_time,     thick_disk_sfr_param,
  -     &        halo_age,                    halo_stars_formation_time
  +     &        thick_disk_max_sfr_relative_time,    thick_disk_sfr_param,
  +     &        halo_age,                    halo_stars_formation_time,
  +     &        thick_disk_birth_init_time,  max_age,
  +     &        halo_birth_init_time
         real :: mass_from_Salpeter_IMF, 
        &        get_normalization_const
         double precision :: coordinate_Theta(MAX_STARS_COUNT),
  
positive.code:       mrep = mrep * massReductionFactor
        normal_mrep = mrep
        burst_mrep = mrep * BURST_FORMATION_FACTOR
        burst_init_time = thick_disk_age - burst_age
  
        thin_disk_birth_init_time = thick_disk_age - thin_disk_age
  C     This can be easily proved by taking derivative from
  C     y = t * exp(-t / tau)
        thick_disk_max_sfr_time = thick_disk_sfr_param
        thick_disk_max_sfr = (thick_disk_max_sfr_time 
       &                      * exp(-thick_disk_max_sfr_time 
       &                              / thick_disk_sfr_param))
  
        write(6,*) '      Mass reduction factor = ', massReductionFactor
  
negative.code:      &                total_generated_mass_in_bin + star_mass
                end if
  
  C             TODO: add halo stars here
                if (disk_belonging(stars_count) == 1) then
                    scale_height = THIN_DISK_SCALE_HEIGHT_KPC
                else if (disk_belonging(stars_count) == 2) then
  