anchor.diff:       &                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
  
positive.code:      &        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
        real :: mass_from_Salpeter_IMF, 
       &        get_normalization_const
        double precision :: coordinate_Theta(MAX_STARS_COUNT),
  
negative.code:       mrep = mrep * massReductionFactor
        normal_mrep = mrep
        burst_mrep = mrep * BURST_FORMATION_FACTOR
  
        max_age = max(thin_disk_age, thick_disk_age, halo_age)
        burst_init_time = max_age - burst_age
        thin_disk_birth_init_time = max_age - thin_disk_age
        thick_disk_birth_init_time = max_age - thick_disk_age
        halo_birth_init_time = max_age - halo_age
  C     This can be easily proved by taking derivative from
  C     y = t * exp(-t / tau)
        thick_disk_max_sfr_relative_time = thick_disk_sfr_param
        thick_disk_max_sfr = (thick_disk_max_sfr_relative_time 
       &                      * exp(-thick_disk_max_sfr_relative_time 
       &                              / thick_disk_sfr_param))
  
        write(6,*) '      Mass reduction factor = ', massReductionFactor
  