use strict;

# transmission dynamics of Dengue-like (vector-borne & mutually cross-reactive)
# viruses.



######################################################
# Definitions of global variables                    #
######################################################

our $script_file = "den_decreasing_transmission_102.pl";

# modifications
# from rev. 98 to 99: The vaccine does not confer cross-protection.
# from rev. 99 to 100: 
#                Natural settings were indroduced as in:
#                  Age-specific birthrate and mortality are input from the file.
#                  Separation between regions were represented by $SEPARATION
#                  Some parameters (properi, nec) were generated as probability distribution function
# from rev. 100 to 101:
#              (1) $VAC_CROSS_FLAG was introduced. 
#                  $VAC_CROSS_FLAG=1: the vaccine induces cross-serotype immunity (e.g. WRAIR's classical vaccine).
#                                 =0: 
#              (2) Effective R was calcualted.
#                  Effective R was defined as the number of viremic individuals in a time-step divided by 
#                  the number of viremic individuals in the previous time-step.
# from rev. 101 to 102: if an individual acquires antibodies to L serotypes,
#                   he/she acquires antibodies to other serotypes, too.

######################################################
# files                                              #
######################################################

#our $file_age_specific_pop = "Thailand_Pop_2005.csv";
our $file_age_specific_pop = "Thailand_Pop_1960.csv";

our $file_age_specific_death = "Thailand_death_2005.csv";


######################################################
# define constant                                    #
######################################################

our $DENV1 = 1;
our $DENV2 = 2;
our $DENV3 = 3;
our $DENV4 = 4;

# non_enhancing vs enhancing

our $VACCINE_NO_VACCINE = 0;
our $VACCINE_NON_ENHANCING = 1;
our $VACCINE_ENHANCING = 2;

#
our $DEAD = 1;
our $ALIVE = 2;

our $FEMALE = 0;
our $MALE = 1;

our $pi = atan2(1, 1) * 4; 

#   perfect sterile immunity to homo-serotype, and transient cross-protection to hetero-serotype

# vaccine BRAND

our $BRAND_WRAIR_2DOSE      = 1;  # WRAIR, 2-dose 0-6-month, formulation 17
our $BRAND_WRAIR_2DOSE_60   = 2;  # WRAIR, 2-dose 0-6-month, formulation 17, seroconversion rate is reduced to 60%
our $BRAND_CHIMERIVAX_BETTER_3DOSE = 3; # Chimerivax, 3-dose 0-6-12-month, better seroconversion reported by Morrison

our $BRAND_CHIMERIVAX_WORSE_3DOSE = 4; # Chimerivax, 3-dose 0-6-12-month, worse seroconversion reported by Poo
our $BRAND_CHIMERIVAX_WORSE_2DOSE = 5; # Chimerivax, 2-dose 0-6-month, worse seroconversion
our $BRAND_CHIMERIVAX_WORSE_2DOSE_60 = 6;  # Chimerivax, 2-dose 0-6-month, worse seroconversion, seroconversion rate is reduced to 60%

our $BRAND_CHIMERIVAX_BETTER_3DOSE_80 = 7;  # Chimerivax, 3-dose 80%
our $BRAND_CHIMERIVAX_BETTER_3DOSE_60 = 8;  # Chimerivax, 3-dose 60%
our $BRAND_CHIMERIVAX_BETTER_3DOSE_40 = 9;  # Chimerivax, 3-dose 40%
our $BRAND_CHIMERIVAX_BETTER_3DOSE_20 =10;  # Chimerivax, 3-dose 40%
our $BRAND_CHIMERIVAX_BETTER_3DOSE_0 = 11;  # Chimerivax, 3-dose 0%
our $BRAND_CHIMERIVAX_BETTER_3DOSE_10 = 14;  # Chimerivax, 3-dose 0%

our $BRAND_CHIMERIVAX_BETTER_2DOSE_40 = 12;  # Chimerivax, 3-dose 40%
our $BRAND_CHIMERIVAX_BETTER_2DOSE_20 = 13;  # Chimerivax, 3-dose 40%

#####################################################
#   Fixed parameter                                 #
#####################################################
our $MINIMUM_RISK = 0.0001; # A virus serotype is revived with this probability
                        # when it is extinct.
                        # The value of $MINIMUM_RISK slightly affects the result.
                        # A larger value of $MINIMUM_RISK sometimes reduces
                        # the total incidence.  
                        # However, this does not result in a qualitative difference.

our $SIM_YEAR = 150;

our $VACCINATION_START_YEAR = 100;

our $transient_period = 10;  # years

our $MAXIMUM_AGE_IN_YEAR = 80;		# maximum age to be used for initialization


our $SEROTYPE = 4;

our $NUM_OF_REGIONS = 3;				# each region is separated from other regions, with a factor of $SEPARATION

our $NATURAL_AGE_STRUCTURE_FLAG = 1;

our $POPULATION;

unless ( $NATURAL_AGE_STRUCTURE_FLAG > 0 ){
	$POPULATION = 100000;
}
else{
	$POPULATION = 0;  # Initialize population emulating the natural demographic structure
}

our $current_population_pointer = $POPULATION;

our $SECONDARY_DHF_PROB = 0.04;

our $PRIMARY_DHF_PROB = 0.002;

our $VACCINATION_AGE = 2; # by years
our $VACCINATION_INTERVAL = 0.5; # by years


our $CUTOFF_1 = 5;	# cutoff age between age classes 1 and 2.
our $CUTOFF_2 = 25;	# cutoff age between age classes 2 and 3.

our $AGE_DEPEND = 0; # age dependency in manifesting DHF. 0: no age dependency

###################################################################

our @SEROCON_RATE;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_WRAIR_2DOSE][1] = 0.17;
$SEROCON_RATE[ $DENV2 ][$BRAND_WRAIR_2DOSE][1] = 0.65;
$SEROCON_RATE[ $DENV3 ][$BRAND_WRAIR_2DOSE][1] = 0.3;
$SEROCON_RATE[ $DENV4 ][$BRAND_WRAIR_2DOSE][1] = 0.65;

$SEROCON_RATE[ $DENV1 ][$BRAND_WRAIR_2DOSE][2] = 0.63;
$SEROCON_RATE[ $DENV2 ][$BRAND_WRAIR_2DOSE][2] = 1;
$SEROCON_RATE[ $DENV3 ][$BRAND_WRAIR_2DOSE][2] = 0.73;
$SEROCON_RATE[ $DENV4 ][$BRAND_WRAIR_2DOSE][2] = 0.83;

$SEROCON_RATE[ $DENV1 ][$BRAND_WRAIR_2DOSE][3] = 0;
$SEROCON_RATE[ $DENV2 ][$BRAND_WRAIR_2DOSE][3] = 0;
$SEROCON_RATE[ $DENV3 ][$BRAND_WRAIR_2DOSE][3] = 0;
$SEROCON_RATE[ $DENV4 ][$BRAND_WRAIR_2DOSE][3] = 0;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_WRAIR_2DOSE_60][1] = 0.17 * 0.6;
$SEROCON_RATE[ $DENV2 ][$BRAND_WRAIR_2DOSE_60][1] = 0.65 * 0.6;
$SEROCON_RATE[ $DENV3 ][$BRAND_WRAIR_2DOSE_60][1] = 0.3 * 0.6;
$SEROCON_RATE[ $DENV4 ][$BRAND_WRAIR_2DOSE_60][1] = 0.65 * 0.6;

$SEROCON_RATE[ $DENV1 ][$BRAND_WRAIR_2DOSE_60][2] = 0.63 * 0.6;
$SEROCON_RATE[ $DENV2 ][$BRAND_WRAIR_2DOSE_60][2] = 1 * 0.6;
$SEROCON_RATE[ $DENV3 ][$BRAND_WRAIR_2DOSE_60][2] = 0.73 * 0.6;
$SEROCON_RATE[ $DENV4 ][$BRAND_WRAIR_2DOSE_60][2] = 0.83 * 0.6;

$SEROCON_RATE[ $DENV1 ][$BRAND_WRAIR_2DOSE_60][3] = 0;
$SEROCON_RATE[ $DENV2 ][$BRAND_WRAIR_2DOSE_60][3] = 0;
$SEROCON_RATE[ $DENV3 ][$BRAND_WRAIR_2DOSE_60][3] = 0;
$SEROCON_RATE[ $DENV4 ][$BRAND_WRAIR_2DOSE_60][3] = 0;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][1] = 0.3;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][1] = 0.73;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][1] = 0.33;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][1] = 0.88;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][2] = 0.89;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][2] = 0.96;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][2] = 0.99;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][2] = 0.92;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][3] = 1;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][3] = 1;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][3] = 1;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE][3] = 1;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][1] = 0.15;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][1] = 0.32;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][1] = 0.50;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][1] = 0.67;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][2] = 0.41;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][2] = 0.63;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][2] = 0.44;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][2] = 0.61;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][3] = 0.54;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][3] = 0.40;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][3] = 0.64;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_3DOSE][3] = 0.15;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][1] = 0.15;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][1] = 0.32;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][1] = 0.5;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][1] = 0.67;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][2] = 0.41;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][2] = 0.63;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][2] = 0.44;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][2] = 0.61;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][3] = 0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][3] = 0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][3] = 0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_2DOSE][3] = 0;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][1] = 0.09;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][1] = 0.19;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][1] = 0.30;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][1] = 0.40;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][2] = 0.25;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][2] = 0.38;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][2] = 0.26;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][2] = 0.37;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][3] = 0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][3] = 0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][3] = 0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_WORSE_2DOSE_60][3] = 0;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][1] = 0.3*0.8;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][1] = 0.73*0.8;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][1] = 0.33*0.8;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][1] = 0.88*0.8;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][2] = 0.89*0.8;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][2] = 0.96*0.8;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][2] = 0.99*0.8;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][2] = 0.92*0.8;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][3] = 1*0.8;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][3] = 1*0.8;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][3] = 1*0.8;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_80][3] = 1*0.8;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][1] = 0.3*0.6;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][1] = 0.73*0.6;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][1] = 0.33*0.6;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][1] = 0.88*0.6;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][2] = 0.89*0.6;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][2] = 0.96*0.6;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][2] = 0.99*0.6;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][2] = 0.92*0.6;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][3] = 1*0.6;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][3] = 1*0.6;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][3] = 1*0.6;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_60][3] = 1*0.6;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][1] = 0.3*0.4;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][1] = 0.73*0.4;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][1] = 0.33*0.4;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][1] = 0.88*0.4;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][2] = 0.89*0.4;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][2] = 0.96*0.4;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][2] = 0.99*0.4;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][2] = 0.92*0.4;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][3] = 1*0.4;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][3] = 1*0.4;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][3] = 1*0.4;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_40][3] = 1*0.4;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][1] = 0.3*0.2;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][1] = 0.73*0.2;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][1] = 0.33*0.2;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][1] = 0.88*0.2;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][2] = 0.89*0.2;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][2] = 0.96*0.2;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][2] = 0.99*0.2;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][2] = 0.92*0.2;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][3] = 1*0.2;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][3] = 1*0.2;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][3] = 1*0.2;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_20][3] = 1*0.2;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][1] = 0.3*0.1;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][1] = 0.73*0.1;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][1] = 0.33*0.1;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][1] = 0.88*0.1;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][2] = 0.89*0.1;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][2] = 0.96*0.1;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][2] = 0.99*0.1;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][2] = 0.92*0.1;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][3] = 1*0.1;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][3] = 1*0.1;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][3] = 1*0.1;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_10][3] = 1*0.1;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][1] = 0.3*0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][1] = 0.73*0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][1] = 0.33*0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][1] = 0.88*0;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][2] = 0.89*0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][2] = 0.96*0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][2] = 0.99*0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][2] = 0.92*0;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][3] = 1*0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][3] = 1*0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][3] = 1*0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_3DOSE_0][3] = 1*0;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][1] = 0.3*0.4;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][1] = 0.73*0.4;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][1] = 0.33*0.4;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][1] = 0.88*0.4;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][2] = 0.89*0.4;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][2] = 0.96*0.4;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][2] = 0.99*0.4;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][2] = 0.92*0.4;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][3] = 0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][3] = 0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][3] = 0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_40][3] = 0;

###################################################################

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][1] = 0.3*0.2;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][1] = 0.73*0.2;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][1] = 0.33*0.2;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][1] = 0.88*0.2;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][2] = 0.89*0.2;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][2] = 0.96*0.2;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][2] = 0.99*0.2;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][2] = 0.92*0.2;

$SEROCON_RATE[ $DENV1 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][3] = 0;
$SEROCON_RATE[ $DENV2 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][3] = 0;
$SEROCON_RATE[ $DENV3 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][3] = 0;
$SEROCON_RATE[ $DENV4 ][$BRAND_CHIMERIVAX_BETTER_2DOSE_20][3] = 0;
###################################################################
our @NUM_OF_DOSES;

$NUM_OF_DOSES[ $BRAND_WRAIR_2DOSE ]     = 2;  # WRAIR, 2-dose 0-6-month, formulation 17
$NUM_OF_DOSES[ $BRAND_WRAIR_2DOSE_60 ]  = 2;  # WRAIR, 2-dose 0-6-month, formulation 17, seroconversion rate is reduced to 60%
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_3DOSE ] = 3; # Chimerivax, 3-dose 0-6-12-month, better seroconversion reported by Morrison

$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_WORSE_3DOSE ]= 3; # Chimerivax, 3-dose 0-6-12-month, worse seroconversion reported by Poo
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_WORSE_2DOSE ]= 2; # Chimerivax, 2-dose 0-6-month, worse seroconversion
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_WORSE_2DOSE_60 ] = 2;  # Chimerivax, 2-dose 0-6-month, worse seroconversion, seroconversion rate is reduced to 60%

$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_3DOSE_80 ]= 3;  # Chimerivax, 3-dose 80%
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_3DOSE_60 ]= 3;  # Chimerivax, 3-dose 60%
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_3DOSE_40 ]= 3;  # Chimerivax, 3-dose 40%
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_3DOSE_20 ]=3;  # Chimerivax, 3-dose 40%
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_3DOSE_0 ]= 3;  # Chimerivax, 3-dose 0%

$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_2DOSE_40 ]= 2;  # Chimerivax, 3-dose 40%
$NUM_OF_DOSES[ $BRAND_CHIMERIVAX_BETTER_2DOSE_20 ]= 2;  # Chimerivax, 3-dose 40%

###################################################################

####################################################################
# parameters given as ARGVs                                        #
####################################################################

our @paraname;
our $COVERAGE;
our $BRAND;

our $MASSVAX = 0;

our $initial_R0;
our $reduced_R0;
our $PROTECTIVE_PERIOD;			# by half months
our $NECESSARY_SEROTYPES;
our $ENHANCED_VIREMIA;
our $FIXED_TFR;
our $SEPARATION;		# 1: perfect separation, 0: no separation between regions.
our $VACCINATION_TYPE;
our $MINIMUM_TRANSMISSION = 0.8;		# represents seasonality.  If $MINIMUM_TRANSMISSION=0, no transmission during winter 

our $PDF_flag;		 	# if this flag is on, probability distribution function will be generated 
										# with mean=ARGV, and SD=1/mean, for the following variables:
										# $PROTECTIVE_PERIOD;			# by half months
										# $ENHANCED_VIREMIA;
our $VAC_CROSS_flag; # if this flag is on, the vaccine confers cross-serotype immunity.

###########################################################################
our $monthly_file_name;
our $annual_file_name;
our $setting_file;
our $age_specific_report_file;

our @R0_of_year;

##################################################################################################
our @indiv_vaccinated_count;	# how many times a vaccine was vaccinated.
our @indiv_vaccinated_timestep;	# the time step at which the individual was vaccinated.

our @indiv_protective_sero_specific_immunity_flag;	# the time step at which the individual acquired serotype-specific Ab.
our @indiv_protective_epitope_sero_specific_inoculation_time_step;	

our @indiv_enhancement_by_wt_flag;
our @indiv_enhancement_by_vaccine_flag;
our @indiv_age_by_half_month;	# age counted by the number of time steps.
our @indiv_gender;

our @indiv_region; # geographical region [1-3]
our @indiv_alive_or_dead; 
##################################################################################################

our @num_of_viremic_individuals_in_timestep;	# $num_of_viremic_individuals_in_timestep[$hmonth]

our @annual_vaccination_per_dose;	# number of vaccinations per dose.
our @annual_vaccination_total;		# total number of vaccinations.

our @monthly_denominator;
our @vaccinated_denominator;	# number of individuals who have been vaccinated at least once.
our @non_vaccinated_denominator;	# number of individuals who have never been vaccinated.

our @monthly_DHF_incidence;
our @monthly_viremic;
our @monthly_inoculation_rate;
our @monthly_total_DHF_incidence;

our @percent_DHF;		# percentage of DHF over total number of inoculations at a month.

our @monthly_primary_DHF_incidence;
our @monthly_DF_incidence;

our @monthly_DHF_incidence_in_non_vaccinated;
our @monthly_DHF_incidence_in_vaccinated;

our @monthly_DF_incidence_in_non_vaccinated;
our @monthly_DF_incidence_in_vaccinated;

our @monthly_DHF_inc_in_non_vac_per_serotype;
our @monthly_DHF_inc_in_vac_per_serotype;

our @risk;	# force of infection (at each time step)
our @extinct_cnt;
our @monthly_risk;
our @weighted_viremia;	# number of viremic individuals for each serotype at a time-step, where DHF enhances viremia.
our @total_inoculation_number;
our @number_of_viremic_individuals;
our @monthly_DF_number;
our @monthly_DF_number_in_vaccinated;
our @monthly_DF_number_in_non_vaccinated;
our @annual_sum_of_age_of_DHF;
our @DHF_number;					# for each serotype, monthly
our @total_DHF_number;		# for all serotypes, monthly
our @primary_DHF_number;	# for all serotypes, monthly
our @age_specific_DHF_number;	# annual number of DHF cases for each age-class

our @DHF_number_in_vaccinated;								# monthly
our @DHF_number_in_vaccinated_serotype;				# monthly
our @DHF_number_in_non_vaccinated;						# monthly
our @DHF_number_in_non_vaccinated_serotype;		# monthly

our @annual_DHF_number_due_to_vaccine;
our @annual_DHF_number_due_to_vaccine_only;

our @annual_inoculation_rate;
our @annual_total_DHF_number;
	
our @annual_DHF_incidence_in_non_vaccinated;
our @annual_DF_incidence_in_non_vaccinated;
our @annual_DHF_incidence_in_vaccinated;
our @annual_DF_incidence_in_vaccinated;
our @annual_total_DHF_incidence;
our @annual_DF_incidence;
our @annual_primary_DHF_incidence;
our @annual_DHF_inc_in_non_vac_per_serotype;
our @annual_DHF_inc_in_vac_per_serotype;
our @annual_risk;			# force of infection
our @annual_DHF_incidence;
our @annual_viremic;
our @annual_total_population;
our @annual_pop_region;		

our @annual_sum_of_effective_R;
our @annual_cnt_of_effective_R;


our @age_specific_DHF_incidence;


######################################################
# demographic 

our $fixed_TFR;

our $output_TFR;

our $population_divide_factor = 1000; # divide the actual population with this.

our @sim_age_specific_pop; # age_specific popoulation used in the simulation

our @actual_age_specific_pop; # actual population: $sim_age_specific_pop/$population_divide_factor

our @age_specific_death; 		# $age_specific_death[ $age_class_cnt ]

our @age_specfic_annual_mortality_rate;


our @age_specific_birthrate;

our $max_age_class = 15;

our $min_fertile_age_class = 4;		# 15-19
our $max_fertile_age_class = 10;	# 45-49

our @tfr_coeff;				# Age-class-specific birthrate will be reconstructed from the total fertility rate (TFR)
our @tfr2_coeff;			# by using second-order regression equations.  
our @cons;						#
#######################################################
# regression coefficients to predict age-class-specific birthrate
$tfr_coeff[4] = 11.46051;
$tfr2_coeff[4] = 1.601075;
$cons[4] = 3.536474;
#######################################################
$tfr_coeff[5] = 50.81174;
$tfr2_coeff[5] = -1.260105;
$cons[5] = 4.477458;
#######################################################
$tfr_coeff[6] = 46.98442;
$tfr2_coeff[6] = -.998202;
$cons[6] = 29.07806;
#######################################################
$tfr_coeff[7] = 39.16125;
$tfr2_coeff[7] = -.3205762;
$cons[7] = 11.31498;
#######################################################
$tfr_coeff[8] = 33.48249;
$tfr2_coeff[8] = -.1705431;
$cons[8] = -22.31366;
#######################################################
$tfr_coeff[9] = 16.14774;
$tfr2_coeff[9] = .340163;
$cons[9] = -20.81537;
#######################################################
$tfr_coeff[10] = 1.951835;
$tfr2_coeff[10] = .8081902;
$cons[10] = -5.277931;
#######################################################


&main;
######################################################
# main                                               #
######################################################
sub main {

	&Request_parameters;

	&Save_settings;

	&Open_outfiles;

	&SetAnnualR0;

	&SetAge;

	&Age_specific_report;

	print "simulation is started\n";

	for (my $hmonth_cnt=0; $hmonth_cnt < $SIM_YEAR * 24; $hmonth_cnt++ ){
		print "\n";
		print "$hmonth_cnt";
		
		my $year = int ($hmonth_cnt / 24);

		if ( $hmonth_cnt % 2 == 0 ){ # recorded once per month
			&Record_immune_state( $hmonth_cnt );
		}
		
		if ( $hmonth_cnt % 24 == 0 ){
			&Record_annually( $hmonth_cnt );
		}
		

		&Increment_age ($hmonth_cnt);

		&Transmit_virus_to_environment ($hmonth_cnt);

		&Transmit_virus_to_humans ($hmonth_cnt);

		&NaturalBirthDeath;
		
		&Calc_Effective_R($hmonth_cnt);
	}

	&PrintMonthlyReport;

	&PrintAnnualReport;

	&Close_outfiles;
}
######################################################
# END OF main                                        #
######################################################


######################################################
# Request_parameters                                 #
######################################################
sub Request_parameters {


	my $trial_sequence = $ARGV[0];
	$paraname[0] = "trial_sequence";
	print "para0:trial_sequence :$trial_sequence\n";
	if ( $trial_sequence eq "" ){
		print "enter arbitrary sequence to identify this simulation.";
		exit;
	}

	$initial_R0 = $ARGV[1];
	$paraname[1] = "initial_R0";
	print "para1: initial_R0 [>0] :$initial_R0\n";
	if ( $initial_R0 eq "" ){
		print "enter initial_R0!";
		exit;
	}

	$reduced_R0 = $ARGV[2];
	$paraname[2] = "reduced_R0";
	print "para2: reduced_R0 [>0] :$reduced_R0\n";
	if ( $reduced_R0 eq "" ){
		print "enter reduced_R0!";
		exit;
	}

	$PROTECTIVE_PERIOD = $ARGV[3];
	$paraname[3] = "PROTECTIVE_PERIOD";
	print "para3: PROTECTIVE_PERIOD: $PROTECTIVE_PERIOD\n";
	if ( $PROTECTIVE_PERIOD eq "" ){
		print "enter PROTECTIVE_PERIOD by half month (i.e. 24 for one year)!\n";
		exit;
	}

	$NECESSARY_SEROTYPES = $ARGV[4];
	$paraname[4] = "NECESSARY_SEROTYPES";
	print "para4: NECESSARY_SEROTYPES: $NECESSARY_SEROTYPES\n";
	if ( $NECESSARY_SEROTYPES eq "" ){
		print "enter NECESSARY_SEROTYPES (number of serotypes which confer total resistance to DHF from the remaining serotypes)!\n";
		exit;
	}

	$ENHANCED_VIREMIA = $ARGV[5];
	$paraname[5] = "ENHANCED_VIREMIA";
	print "para5: ENHANCED_VIREMIA: $ENHANCED_VIREMIA\n";
	if ( $ENHANCED_VIREMIA eq "" ){
		print "enter ENHANCED_VIREMIA [>=1] (1: no enhancement in transmissibility during DHF)!\n";
		exit;
	}

	$FIXED_TFR = $ARGV[6];
	$paraname[6] = "FIXED_TFR";
	print "para6: FIXED_TFR: $FIXED_TFR\n";
	if ( $FIXED_TFR eq "" ){
		print "enter FIXED TFR\n";
		exit;
	}

	$SEPARATION = $ARGV[7];
	$paraname[7] = "SEPARATION";
	print "para7: SEPARATION: $SEPARATION\n";
	if ( $SEPARATION eq "" ){
		print "enter SEPARATION\n";
		exit;
	}

	$MINIMUM_TRANSMISSION = $ARGV[8];
	$paraname[8] = "MINIMUM_TRANSMISSION";
	print "para8: MINIMUM_TRANSMISSION: $MINIMUM_TRANSMISSION\n";
	if ( $MINIMUM_TRANSMISSION eq "") {
		print "enter MINIMUM_TRANSMISSION\n";
		exit;
	}

	$VACCINATION_TYPE = $ARGV[9];
	$paraname[9] = "VACCINATION_TYPE";
	print "para9: VACCINATION_TYPE: $VACCINATION_TYPE\n";
	if ( $VACCINATION_TYPE eq "" ){
		print "enter VACCINATION_TYPE [0, 1, 2] (0: no vaccination, 1: non-enhancing, 2: enhancing vaccine\n";
		exit;
	}

	$COVERAGE = $ARGV[10];
	$paraname[10] = "COVERAGE";
	print "para10: COVERAGE: $COVERAGE\n";
	if ( $COVERAGE eq "" ){
		print "enter COVERAGE [0--1]\n";
		exit;
	}
	if ( $COVERAGE > 1){
		print "COVERAGE should be [0--1]\n";
		exit;
	}

	$BRAND = $ARGV[11];
	$paraname[11] = "BRAND";
	print "para11: BRAND: $BRAND\n";
	if ( $BRAND eq "" ){
		print "enter BRAND of the vaccine regimen \n";
		exit;
	}

	$VAC_CROSS_flag = $ARGV[12];
	$paraname[12] = "VAC_CROSS_flag";
	print "para12: VAC_CROSS_flag: $VAC_CROSS_flag\n";
	if ( $VAC_CROSS_flag eq "" ){
		print "enter VAC_CROSS_flag\n";
		exit;
	}

	$PDF_flag = $ARGV[13];
	$paraname[13] = "PDF_flag";
	print "para13: PDF_flag: $PDF_flag\n";
	if ( $BRAND eq "" ){
		print "enter PDF_flag\n";
		exit;
	}
	

	my $file_suffix = "sq_${trial_sequence}_iR0_${initial_R0}_rR0_${reduced_R0}_pp_${PROTECTIVE_PERIOD}_L_${NECESSARY_SEROTYPES}_enh_${ENHANCED_VIREMIA}_tf_${FIXED_TFR}_sp_${SEPARATION}_mt_${MINIMUM_TRANSMISSION}_vt_${VACCINATION_TYPE}_cv_${COVERAGE}_br_${BRAND}_vcros_${VAC_CROSS_flag}_pdf_${PDF_flag}";

	$monthly_file_name = "Simulation_Results\\DLV_Monthly_" . $file_suffix . ".csv";

	$annual_file_name = "Simulation_Results\\DLV_Annual_" . $file_suffix . ".csv";

	$setting_file = "Simulation_Results\\DLV_Settings_for_" . $file_suffix . ".txt";

	$age_specific_report_file = "Simulation_Results\\DLV_Age_report_for_" . $file_suffix . ".csv";


}
######################################################
# END OF Request_parameters                          #
######################################################

######################################################
# Save_settings                                      #
######################################################
sub Save_settings {
	if (open (TEMP2, $setting_file) ){
		print "Overwriting a existing file!!\n";
		close (TEMP2);
		exit;
	}

	if ( open (SETTING,">$setting_file" ) ){
		print SETTING "*************************************************\n";
		for (my $argcnt = 0; $argcnt < @ARGV; $argcnt++){
			print SETTING "ARGV[$argcnt]:$paraname[$argcnt] = $ARGV[$argcnt]\n";
		}
		print SETTING "*************************************************\n";

		if (open (PRG, "$script_file") ){
			while (<PRG>){
				print SETTING;
			}
			close (PRG);
		}
		else {
			print "programme file not opened\n";
			exit;
		}
		close (SETTING);
	}
	else {
		print "setting file not opened\n";
		exit;
	}
}
######################################################
# END OF Save_settings                               #
######################################################

######################################################
# Open_outfiles                                      #
######################################################
sub Open_outfiles{

	unless (open MONTHLY, ">$monthly_file_name" ){
		print "MONTHLY not opened\n";
		exit;
	}

	unless (open ANNUAL, ">$annual_file_name" ){
		print "ANNUAL not opened\n";
		exit;
	}
}
######################################################
# END OF Open_outfiles                               #
######################################################

######################################################
# Close_outfiles                                     #
######################################################
sub Close_outfiles{
	close (MONTHLY);
	close (ANNUAL);
}
######################################################
# END OF Close_outfiles                              #
######################################################

######################################################
# SetAnnualR0                                        #
######################################################
sub SetAnnualR0 {
	my $year;

	my $Vector_control_start_year = $VACCINATION_START_YEAR;


	for ( my $year = 0; $year < $Vector_control_start_year; $year++ ){
		$R0_of_year[ $year ] = $initial_R0;
	}
	for ( $year = $Vector_control_start_year; $year < ($Vector_control_start_year + $transient_period); $year++ ){
		$R0_of_year[ $year ] = ($initial_R0*($Vector_control_start_year+$transient_period-$year) + $reduced_R0*($year-$Vector_control_start_year))/$transient_period;
	}
	for ( $year = $Vector_control_start_year + $transient_period; $year < $SIM_YEAR; $year++ ){
		$R0_of_year[ $year ] = $reduced_R0;
	}

}
######################################################
# END OF SetAnnualR0                                 #
######################################################

######################################################
# SetAge                                             #
######################################################
sub SetAge {
	# set initial age

	my $human_cnt;
	my $sero_cnt;

	if ( $NATURAL_AGE_STRUCTURE_FLAG > 0 ){

		&Read_age_pop_file;
		
		&Read_age_specific_death;

		&Calc_age_specific_mortality;		

		&Calc_age_specific_birthrate;

		&Calc_TFR;
		
		$current_population_pointer = &Initialize_population;
		
		print "current_population_pointer: $current_population_pointer\n";

	
	}
	else {	# Flat age structure
		for ($human_cnt=0; $human_cnt < $current_population_pointer; $human_cnt++ ){

			$indiv_vaccinated_count[ $human_cnt ] = 0;
			$indiv_vaccinated_timestep[ $human_cnt ] = "";

			for ( $sero_cnt = 1; $sero_cnt <= $SEROTYPE; $sero_cnt++ ){

				# protective epitope
				$indiv_protective_sero_specific_immunity_flag[ $sero_cnt ][ $human_cnt] = 0;
				$indiv_protective_epitope_sero_specific_inoculation_time_step[ $sero_cnt ][ $human_cnt] = 24 * ($SIM_YEAR+1000); # impossibly large number

			}

			$indiv_enhancement_by_wt_flag[ $human_cnt ] = 0;

			$indiv_enhancement_by_vaccine_flag[ $human_cnt ] = 0;

			$indiv_age_by_half_month[$human_cnt] = ($MAXIMUM_AGE_IN_YEAR * 24) * rand();

		}
	}
}
######################################################
# END OF SetAge                                      #
######################################################

######################################################
# Age_specific_report                                #
######################################################
sub Age_specific_report {
	if (open (AGE_REP, ">$age_specific_report_file" ) ){

		print AGE_REP "class,actual_pop,sim_pop,death,annual_mortality_rate,birthrate\n";

		for (my $age_class_cnt = 1; $age_class_cnt <= $max_age_class; $age_class_cnt++ ){
			print AGE_REP "$age_class_cnt,$actual_age_specific_pop[ $age_class_cnt ],$sim_age_specific_pop[ $age_class_cnt ],$age_specific_death[ $age_class_cnt ],$age_specfic_annual_mortality_rate[ $age_class_cnt ],$age_specific_birthrate[ $age_class_cnt ]\n";
		}
		close(AGE_REP);
	}
}
######################################################
# END OF Age_specific_report                         #
######################################################

	
######################################################
# Read_age_pop_file                                  #
######################################################
sub Read_age_pop_file {

	if (open (AGE_POP, $file_age_specific_pop) ){
		while(<AGE_POP>){
			s/ //g;
			s/\n//;
			if (/^\d+/){
				my @elements = split(/,/, $_);
				my $total_pop = $elements[ 0 ];
				my $pop_temp = 0;
				for ( my $age_class_cnt = 1; $age_class_cnt < $max_age_class; $age_class_cnt++ ){
					#pp0_4	pp5_9	pp10_14	pp15_19	pp20_24	pp25_29	pp30_34	pp35_39	pp40_44	pp45_49	pp50_54	pp55_59	pp60_64	pp65_69
					#1			2			3				4				5				6				7				8				9				10			11			12			13			14
					$actual_age_specific_pop[ $age_class_cnt ] = $elements[ $age_class_cnt ];
					$pop_temp += $elements[ $age_class_cnt ];
				}
				$actual_age_specific_pop[ $max_age_class ] = $total_pop - $pop_temp;  # class 15: 70 years and above
			}
		}
		close(AGE_POP);
		
		for (my $age_class_cnt = 1; $age_class_cnt <= $max_age_class; $age_class_cnt++ ){
			$sim_age_specific_pop[ $age_class_cnt ] = int( $actual_age_specific_pop[ $age_class_cnt ]/$population_divide_factor );
		}
		
	}
	else {
		print "AGE_POP not opened\n";
		exit;
	}
}
######################################################
# END OF Read_age_pop_file                           #
######################################################

######################################################
# Read_age_specific_death                            #
######################################################
sub Read_age_specific_death {
	if (open( AGE_DEATH, $file_age_specific_death) ){
		my $age_class_cnt = 1;
		while(<AGE_DEATH>){
			s/ //g;
			s/\n//;
			if (/^death/){
				my @elements = split(/,/, $_);
				$age_specific_death[ $age_class_cnt ] = $elements[ 1 ];
				$age_class_cnt++;
			}
		}
		close(AGE_DEATH);
	}
	else {
		print "AGE_DEATH not opened\n";
		exit;
	}
}
######################################################
# END OF Read_age_specific_death                     #
######################################################

######################################################
# Calc_age_specific_mortality                        #
######################################################
sub Calc_age_specific_mortality {
	for (my $age_class_cnt = 1; $age_class_cnt <= $max_age_class; $age_class_cnt++ ){
		$age_specfic_annual_mortality_rate[ $age_class_cnt ] = $age_specific_death[ $age_class_cnt ]/$actual_age_specific_pop[ $age_class_cnt ];
	}
}
######################################################
# END OF Calc_age_specific_mortality                 #
######################################################

######################################################
# Calc_age_specific_birthrate                        #
######################################################
sub Calc_age_specific_birthrate {

	for (my $age_class_cnt = $min_fertile_age_class; $age_class_cnt <= $max_fertile_age_class; $age_class_cnt++ ){
		my $birth_per_1000 = $tfr_coeff[ $age_class_cnt ]*$FIXED_TFR + $tfr2_coeff[ $age_class_cnt ]*($FIXED_TFR**2) + $cons[$age_class_cnt];
		$age_specific_birthrate[ $age_class_cnt ] = $birth_per_1000 / 1000;
	}
}
######################################################
# END OF Calc_age_specific_birthrate                 #
######################################################

######################################################
# Calc_TFR                                           #
######################################################
sub Calc_TFR {

	# Reconstruct total fertility rate (TFR) from age-class-specific birthrates.
	
	for (my $age_class_cnt = $min_fertile_age_class; $age_class_cnt <= $max_fertile_age_class; $age_class_cnt++ ){
		$output_TFR += $age_specific_birthrate[ $age_class_cnt ] * 5;
	}
	
	print "output_TFR: $output_TFR\n";
}
######################################################
# END OF Calc_TFR                                    #
######################################################


######################################################
# Initialize_population                              #
######################################################
sub Initialize_population {
	my $human_cnt = 0;
	my $sero_cnt;
	
	for (my $age_class_cnt = 1; $age_class_cnt <= $max_age_class; $age_class_cnt++ ){
		my $min_age_within_class = ($age_class_cnt - 1) * 5;
		my $max_age_within_class = $min_age_within_class + 4;
		my $initial_number_per_age = $sim_age_specific_pop[ $age_class_cnt ]/5 ;

		for (my $age_by_year = $min_age_within_class; $age_by_year <= $max_age_within_class; $age_by_year++ ){
			for (my $temp_counter = 0; $temp_counter < $initial_number_per_age; $temp_counter++ ){

				$indiv_age_by_half_month[$human_cnt] = $age_by_year * 24;  # 24 half months per year

				$indiv_vaccinated_count[ $human_cnt ] = 0;
				
				$indiv_vaccinated_timestep[ $human_cnt ] = "";

				$indiv_gender[$human_cnt] = int(rand(2));		# male or female

				$indiv_enhancement_by_wt_flag[ $human_cnt ] = 0;

				#$indiv_enhancement_by_wt_flag[ $human_cnt ] = int( rand(2) );	# At the initial year, half of the population
																																			# are assumed to be seropositive

				$indiv_enhancement_by_vaccine_flag[ $human_cnt ] = 0;

				$indiv_alive_or_dead[ $human_cnt ] = $ALIVE;

				$indiv_region[ $human_cnt ] = int( rand($NUM_OF_REGIONS) );

				for ( $sero_cnt = 1; $sero_cnt <= $SEROTYPE; $sero_cnt++ ){

					# protective epitope
					$indiv_protective_sero_specific_immunity_flag[ $sero_cnt ][ $human_cnt] = 0;
					
					#$indiv_protective_sero_specific_immunity_flag[ $sero_cnt ][ $human_cnt] = int( rand(2) ); # At the initial year,
																																																		# half of the population
																																																		# are assumed to be seropositive.

					$indiv_protective_epitope_sero_specific_inoculation_time_step[ $sero_cnt ][ $human_cnt] = 24 * ($SIM_YEAR+1000); # impossibly large number
				}

				$human_cnt++;
			}
		}
	}

	return ( $human_cnt );

}
#####################################################################
# END OF Initialize_population                                      #
#####################################################################


######################################################
# Increment_age                                      #
######################################################
sub Increment_age {

	my $hmonth = $_[0];
	my $person;

	for ($person=0; $person<$current_population_pointer; $person++){

		if ( $indiv_alive_or_dead[$person] == $ALIVE ){

			$indiv_age_by_half_month[$person]++;

			if ( $VACCINATION_TYPE > 0){

				if ( $hmonth == $VACCINATION_START_YEAR * 24 ){ # vaccinate start year
					if ( $MASSVAX > 0 ){
						if ( $indiv_age_by_half_month[$person] >= $VACCINATION_AGE * 24  ){
							if ( rand() < $COVERAGE ) {
								$indiv_vaccinated_count[ $person ] = 1;
								&Vaccinate($person, $hmonth, 1);	# first dose
								$indiv_vaccinated_timestep[ $person ] = $hmonth;
							}
						}
					}
				}
				elsif ( $hmonth > $VACCINATION_START_YEAR * 24){ # after the first year of vaccination

					if ( $indiv_vaccinated_count[$person] == 0){
						if ($indiv_age_by_half_month[$person] == $VACCINATION_AGE * 24 ){ 
							if ( rand() < $COVERAGE ){
								$indiv_vaccinated_count[ $person ] = 1;
								&Vaccinate($person, $hmonth, 1); # first dose
								$indiv_vaccinated_timestep[ $person ] = $hmonth;
							}
						}
					}
					elsif ( $indiv_vaccinated_count[$person] == 1 ){ # already vaccinated once
						if ( ($hmonth - $indiv_vaccinated_timestep[ $person ]) == $VACCINATION_INTERVAL * 24 ){ 
							if ( rand() < $COVERAGE ){
								$indiv_vaccinated_count[ $person ] = 2;
								&Vaccinate($person, $hmonth, 2); # second dose
								$indiv_vaccinated_timestep[ $person ] = $hmonth;
							}
						}
					}
					elsif ( $indiv_vaccinated_count[$person] == 2 ){ # already vaccinated twice
						if ( $NUM_OF_DOSES[ $BRAND ] == 3 ){
							if ( ($hmonth - $indiv_vaccinated_timestep[ $person ]) == $VACCINATION_INTERVAL * 24 ){ 
								if ( rand() < $COVERAGE ){
									$indiv_vaccinated_count[ $person ] = 3;
									&Vaccinate($person, $hmonth, 3); # third dose
									$indiv_vaccinated_timestep[ $person ] = $hmonth;
								}
							}
						}
					}
				}
			}
		}
	}
}
######################################################
# End of Increment_age                               #
######################################################

######################################################
# Vaccinate                                          #
######################################################
sub Vaccinate {

	my $person = $_[0];
	my $hmonth = $_[1];
	my $dose = $_[2]; # 1: first dose, 2: second dose, 3: third dose
	my $serocnt;
	my $year = int($hmonth/24);

	if ( $VACCINATION_TYPE == $VACCINE_ENHANCING ){ # protective epitope only, tetra valent

		$indiv_enhancement_by_vaccine_flag[ $person ] = 1;

	}

	$annual_vaccination_per_dose[$dose][$year]++;
	$annual_vaccination_total[$year]++;

	for ( $serocnt = 1; $serocnt <= $SEROTYPE; $serocnt++ ){
		if ( rand() < $SEROCON_RATE[ $serocnt ][ $BRAND ][$dose] ){
			$indiv_protective_sero_specific_immunity_flag[ $serocnt ][ $person ] = 1;

			if ( $VAC_CROSS_flag > 0 ){
				# the vaccine confers cross-serotype protection.
			
				$indiv_protective_epitope_sero_specific_inoculation_time_step[ $serocnt ][ $person ] = $hmonth;
			}

		}
	}
}
######################################################
# END OF Vaccinate                                   #
######################################################

######################################################
# Record_immune_state                                #
######################################################
sub Record_immune_state {

	my $half_month = $_[0];
	my $month = int($half_month/2);
	my $year = int($half_month/24);

	my $person;
	my $age_year;
	my $int_age;

	for ( $person = 0; $person < $current_population_pointer; $person++ ){

		if ( $indiv_alive_or_dead[ $person ] == $ALIVE ){

			$age_year = $indiv_age_by_half_month[ $person ] /24;
			$int_age = int ($age_year);
			#$age_specific_denominator[ $month ][$int_age]++;
			$monthly_denominator[$month]++;

		
			if ( $indiv_vaccinated_count[ $person ] > 0 ){
				$vaccinated_denominator[$month]++;
			}
			else {
				$non_vaccinated_denominator[$month]++;
			}
		}
	}
}
######################################################
# END OF Record_immune_state                         #
######################################################

######################################################
# Record_annually                                    #
######################################################
sub Record_annually { 
	my $half_month = $_[0];
	my $year = int($half_month/24);

	if ( $NATURAL_AGE_STRUCTURE_FLAG > 0 ){
		for (my $person = 0; $person < $current_population_pointer; $person++ ){
			if ( $indiv_alive_or_dead[ $person ] == $ALIVE ){
				$annual_total_population[ $year ]++;
			}
		}
	}
	else {
		$annual_total_population[ $year ] = $POPULATION;
	}	
}
######################################################
# END OF Record_annually                             #
######################################################


######################################################
# PrintMonthlyReport                                 #
######################################################
sub PrintMonthlyReport {

	my $outrec;
	my $month;
	my $year;
	my $R0;
	my $sero;

	$outrec =  "month,R0_of_year,monthly_inoculation_rate,percent_DHF,";
	$outrec .= "risk_D1,risk_D2,risk_D3,risk_D4,";
	$outrec .= "total DHF inc,DHF inc in primary infections,";
	$outrec .= "DF inc,";
	$outrec .= "DHF inc in vaccinated,DHF inc in non_vaccinated,";
	$outrec .= "DHF_D1,DHF_D2,DHF_D3,DHF_D4,";
	$outrec .= "viremic_D1,viremic_D2,viremic_D3,viremic_D4\n";

	print MONTHLY "$outrec";


	for ($month = 0; $month < $SIM_YEAR*12; $month++){
		if ( $monthly_denominator[$month] > 0 ){
			$year = int ( $month / 12 );
			$R0 = $R0_of_year[ $year ];

			for ($sero = 1; $sero <=$SEROTYPE; $sero++ ){

				$monthly_DHF_incidence[$sero][$month] = 100000* $DHF_number[$month][$sero] / $annual_total_population[ $year ];

				$monthly_viremic[$sero][$month] = 100000 * $number_of_viremic_individuals[$sero][$month]/$annual_total_population[ $year ];

			}

			$monthly_inoculation_rate[$month] = 100000*$total_inoculation_number[$month]/$annual_total_population[ $year ];

			$monthly_total_DHF_incidence[$month] = 100000 * $total_DHF_number[ $month ]/$annual_total_population[ $year ];

			if ( $monthly_total_DHF_incidence[$month] > 0 ){
				$percent_DHF[$month] = 100 * $monthly_total_DHF_incidence[$month] / $monthly_inoculation_rate[$month];
			}
			else {
				$percent_DHF[$month] = "";
			}

			$monthly_primary_DHF_incidence[$month] = 100000 * $primary_DHF_number[ $month ]/$annual_total_population[ $year ];

			$monthly_DF_incidence[$month] = 100000 * $monthly_DF_number[ $month ]/$annual_total_population[ $year ];

			if ( $non_vaccinated_denominator[$month] > 0 ){ 
				$monthly_DHF_incidence_in_non_vaccinated[$month] = 100000 * $DHF_number_in_non_vaccinated[ $month ]/$non_vaccinated_denominator[$month];
				$monthly_DF_incidence_in_non_vaccinated[$month] = 100000 * $monthly_DF_number_in_non_vaccinated[ $month ]/$non_vaccinated_denominator[$month];

				for ( $sero = 1; $sero<=$SEROTYPE; $sero++ ){
					$monthly_DHF_inc_in_non_vac_per_serotype[$month][$sero] = 100000 * $DHF_number_in_non_vaccinated_serotype[$month][$sero]/$non_vaccinated_denominator[$month];
				}
			}

			if ( $vaccinated_denominator[$month] > 0 ){
				$monthly_DHF_incidence_in_vaccinated[$month] = 100000 * $DHF_number_in_vaccinated[ $month ]/$vaccinated_denominator[$month];
				$monthly_DF_incidence_in_vaccinated[$month] = 100000 * $monthly_DF_number_in_vaccinated[ $month ]/$vaccinated_denominator[$month];

				for ( $sero = 1; $sero<=$SEROTYPE; $sero++ ){
					$monthly_DHF_inc_in_vac_per_serotype[$month][$sero] = 100000 * $DHF_number_in_vaccinated_serotype[$month][$sero]/$vaccinated_denominator[$month];
				}

				#$monthly_DHF_incidence_due_to_vaccine[$month] = 100000 * $DHF_number_due_to_vaccine[$month] / $vaccinated_denominator[$month];

			}

			$outrec =  "$month,$R0,$monthly_inoculation_rate[$month],$percent_DHF[$month],";
			$outrec .= "$monthly_risk[1][$month],$monthly_risk[2][$month],$monthly_risk[3][$month],$monthly_risk[4][$month],";
			$outrec .= "$monthly_total_DHF_incidence[$month],$monthly_primary_DHF_incidence[$month],";
			$outrec .= "$monthly_DF_incidence[$month],";
			$outrec .= "$monthly_DHF_incidence_in_vaccinated[$month],$monthly_DHF_incidence_in_non_vaccinated[$month],";
			$outrec .= "$monthly_DHF_incidence[1][$month],$monthly_DHF_incidence[2][$month],$monthly_DHF_incidence[3][$month],$monthly_DHF_incidence[4][$month],";
			$outrec .= "$monthly_viremic[1][$month],$monthly_viremic[2][$month],$monthly_viremic[3][$month],$monthly_viremic[4][$month]\n";

			print MONTHLY "$outrec";
		}
		else {
			print "No individuals.\n";
			exit;
		}
	}
}
######################################################
# END OF PrintMonthlyReport                          #
######################################################

######################################################
# PrintAnnualReport                                  #
######################################################
sub PrintAnnualReport {

	my $annual_title;
	my $year;
	my $month;
	my $sero_cnt;
	my $recnt;
	my $annual_mean_age;
	my $annual_IMA_DHF;
	my $annual_rec;
	my $average_effective_R;
	my $DHF_due_to_vaccination_over_total_vaccination;
	my $DHF_due_to_vac_only_over_total_vaccination;
	my $DHF_due_to_vaccination_over_first_dose;
	my $DHF_due_to_vac_only_over_first_dose;

	$annual_title = "year,R0,effective_R,population,inoculation_rate [/100 000],";
	$annual_title .= "mean_age,";
	$annual_title .= "total DHF incidence,DHF in primary infections,";
	$annual_title .= "DF incidence,DF inc in vaccinated,DF inc in non_vaccinated,";
	#$annual_title .= "age_specific_DHF_inc_class1,age_specific_DHF_inc_class2,age_specific_DHF_inc_class3,";
	$annual_title .= "DHF incidence in the vaccinated,DHF incidence in the non_vaccinated,";
	$annual_title .= "risk D1,risk D2,risk D3,risk D4,";
	$annual_title .= "DHF D1,DHF D2,DHF D3,DHF D4,";
	$annual_title .= "viremic D1,viremic D2,viremic D3,viremic D4,";
	$annual_title .= "DHF in vac D1,DHF in vac D2,DHF in vac D3,DHF in vac D4,";
	$annual_title .= "DHF in non_vac D1,DHF in non_vac D2,DHF in non_vac D3,DHF in non_vac D4,";
	$annual_title .= "DHF due to vaccine,";
	$annual_title .= "DHF due to vac only,";
	$annual_title .= "vaccination_1st,vaccination_2nd,vaccination_3rd,";
	$annual_title .= "vaccination_all,";
	$annual_title .= "DHF due to vac over total vac [/100 000],DHF due to vac over first dose [/100 000],";
	$annual_title .= "DHF due to vac only over total vac [/100 000],DHF due to vac only over first dose [/100 000],";
	$annual_title .= "extinct_D1,extinct_D2,extinct_D3,extinct_D4\n";

	print ANNUAL "$annual_title";

	for ($year = 0; $year < $SIM_YEAR; $year++ ){
		$recnt = 0;
		for ( $month = 12*$year; $month <= 12*$year+11; $month++ ){
			if ( $monthly_denominator[$month] > 0 ){
				$recnt++;

				$annual_inoculation_rate[ $year ] += $monthly_inoculation_rate[$month];

				$annual_total_DHF_number[ $year ] += $total_DHF_number[ $month ];

				$annual_DHF_incidence_in_non_vaccinated[$year] += $monthly_DHF_incidence_in_non_vaccinated[$month];
				$annual_DF_incidence_in_non_vaccinated[$year] += $monthly_DF_incidence_in_non_vaccinated[$month];

				$annual_DHF_incidence_in_vaccinated[$year] += $monthly_DHF_incidence_in_vaccinated[$month];
				$annual_DF_incidence_in_vaccinated[$year] += $monthly_DF_incidence_in_vaccinated[$month];

				$annual_total_DHF_incidence[ $year ] += $monthly_total_DHF_incidence[$month];

				$annual_DF_incidence[$year] += $monthly_DF_incidence[$month];

				$annual_primary_DHF_incidence[ $year ] += $monthly_primary_DHF_incidence[$month];

				for ( $sero_cnt = 1; $sero_cnt <= $SEROTYPE; $sero_cnt++ ){

					$annual_DHF_inc_in_non_vac_per_serotype[$year][$sero_cnt] += $monthly_DHF_inc_in_non_vac_per_serotype[$month][$sero_cnt];
					$annual_DHF_inc_in_vac_per_serotype[$year][$sero_cnt] += $monthly_DHF_inc_in_vac_per_serotype[$month][$sero_cnt];

					$annual_risk[$sero_cnt][$year] += $monthly_risk[ $sero_cnt ][ $month ];
					$annual_DHF_incidence[$sero_cnt][$year] += $monthly_DHF_incidence[ $sero_cnt ][ $month ];
					$annual_viremic[ $sero_cnt ][$year] += $monthly_viremic[$sero_cnt][$month];

				}
			}
		}

		if ( $recnt > 0 ){
			if ( $annual_total_DHF_number[ $year ] > 0 ){
				$annual_mean_age = $annual_sum_of_age_of_DHF[ $year ] / $annual_total_DHF_number[ $year ];
			}
			else {
				$annual_mean_age = "";
			}

		}
		
		#$age_specific_DHF_incidence[$year][1] = 100000 * $age_specific_DHF_number[$year][1]/($annual_total_population[ $year ]*$CUTOFF_1/$MAXIMUM_AGE_IN_YEAR);
		#$age_specific_DHF_incidence[$year][2] = 100000 * $age_specific_DHF_number[$year][2]/($annual_total_population[ $year ]*($CUTOFF_2-$CUTOFF_1)/$MAXIMUM_AGE_IN_YEAR);
		#$age_specific_DHF_incidence[$year][3] = 100000 * $age_specific_DHF_number[$year][3]/($annual_total_population[ $year ]*($MAXIMUM_AGE_IN_YEAR-$CUTOFF_2)/$MAXIMUM_AGE_IN_YEAR);

		
		if ( $annual_vaccination_total[$year] > 0 ){
			$DHF_due_to_vaccination_over_total_vaccination = 100000 * $annual_DHF_number_due_to_vaccine[$year]/$annual_vaccination_total[$year];
			$DHF_due_to_vac_only_over_total_vaccination = 100000 *$annual_DHF_number_due_to_vaccine_only[ $year ]/$annual_vaccination_total[$year];
		}
		else { 
			$DHF_due_to_vaccination_over_total_vaccination = "";
			$DHF_due_to_vac_only_over_total_vaccination = "";
		}
		
		if ( $annual_vaccination_per_dose[1][$year] > 0 ){
			$DHF_due_to_vaccination_over_first_dose = 100000 * $annual_DHF_number_due_to_vaccine[$year]/$annual_vaccination_per_dose[1][$year];
			$DHF_due_to_vac_only_over_first_dose = 100000 * $annual_DHF_number_due_to_vaccine_only[$year]/$annual_vaccination_per_dose[1][$year];
		}
		else {
			$DHF_due_to_vaccination_over_first_dose = "";
			$DHF_due_to_vac_only_over_first_dose = "";
		}

		if ( $annual_cnt_of_effective_R[ $year ] > 0 ){

			$average_effective_R = $annual_sum_of_effective_R[ $year ] / $annual_cnt_of_effective_R[ $year ];
		}
		else {
			$average_effective_R = "";
		}

		$annual_rec = "$year,$R0_of_year[ $year ],$average_effective_R,$annual_total_population[ $year ],$annual_inoculation_rate[ $year ],";
		$annual_rec .= "$annual_mean_age,";
		$annual_rec .= "$annual_total_DHF_incidence[ $year ],$annual_primary_DHF_incidence[ $year ],";
		$annual_rec .= "$annual_DF_incidence[ $year ],$annual_DF_incidence_in_vaccinated[$year],$annual_DF_incidence_in_non_vaccinated[$year],";
		#$annual_rec .= "$age_specific_DHF_incidence[$year][1],$age_specific_DHF_incidence[$year][2],$age_specific_DHF_incidence[$year][3],";
		$annual_rec .= "$annual_DHF_incidence_in_vaccinated[ $year ],$annual_DHF_incidence_in_non_vaccinated[ $year ],";
		$annual_rec .= "$annual_risk[1][$year],$annual_risk[2][$year],$annual_risk[3][$year],$annual_risk[4][$year],";
		$annual_rec .= "$annual_DHF_incidence[1][$year],$annual_DHF_incidence[2][$year],$annual_DHF_incidence[3][$year],$annual_DHF_incidence[4][$year],";
		$annual_rec .= "$annual_viremic[1][$year],$annual_viremic[2][$year],$annual_viremic[3][$year],$annual_viremic[4][$year],";
		$annual_rec .= "$annual_DHF_inc_in_vac_per_serotype[$year][1],$annual_DHF_inc_in_vac_per_serotype[$year][2],$annual_DHF_inc_in_vac_per_serotype[$year][3],$annual_DHF_inc_in_vac_per_serotype[$year][4],";
		$annual_rec .= "$annual_DHF_inc_in_non_vac_per_serotype[$year][1],$annual_DHF_inc_in_non_vac_per_serotype[$year][2],$annual_DHF_inc_in_non_vac_per_serotype[$year][3],$annual_DHF_inc_in_non_vac_per_serotype[$year][4],";
		$annual_rec .= "$annual_DHF_number_due_to_vaccine[$year],";
		$annual_rec .= "$annual_DHF_number_due_to_vaccine_only[$year],";
		$annual_rec .= "$annual_vaccination_per_dose[1][$year],$annual_vaccination_per_dose[2][$year],$annual_vaccination_per_dose[3][$year],";
		$annual_rec .= "$annual_vaccination_total[$year],";
		$annual_rec .= "$DHF_due_to_vaccination_over_total_vaccination,$DHF_due_to_vaccination_over_first_dose,";
		$annual_rec .= "$DHF_due_to_vac_only_over_total_vaccination,$DHF_due_to_vac_only_over_first_dose,";
		$annual_rec .= "$extinct_cnt[1][$year],$extinct_cnt[2][$year],$extinct_cnt[3][$year],$extinct_cnt[4][$year]";
		print ANNUAL "$annual_rec\n";

	}
}
######################################################
# END OF PrintAnnualReport                           #
######################################################

######################################################
#  Transmit_virus_to_environment                     #
######################################################
sub Transmit_virus_to_environment {

	my $hmon = $_[0];
	my $year = int ($hmon / 24);
	my $month = int ($hmon / 2);
	my $region_pnt;
	
	# minimum: $MINIMUM_TRANSMISSION, mean: 1, maximum: 1+(1-$MINIMUM_TRANSMISSION) = 2-$MINIMUM_TRANSMISSION
	my $seasonality = -cos(2 * $pi * $hmon / 24)*(1-$MINIMUM_TRANSMISSION) + 1;

	# (1 - $MINIMUM_TRANSMISSION) represents the strength of seasonality.
	 
	for (my $sero_cnt=1; $sero_cnt <= $SEROTYPE; $sero_cnt++ ){

		for ($region_pnt = 0; $region_pnt < $NUM_OF_REGIONS; $region_pnt++ ){
			#risk: probability of viral inoculation.

			$risk[$sero_cnt][$hmon][$region_pnt] = $R0_of_year[ $year ] * ($weighted_viremia[$sero_cnt][$region_pnt] / $annual_total_population[ $year ] ) * $seasonality;

			for (my $other_region = 0; $other_region < $NUM_OF_REGIONS; $other_region++ ){

				if ( $other_region != $region_pnt ){
					$risk[$sero_cnt][$hmon][$region_pnt] += (1 - $SEPARATION) * $R0_of_year[ $year ] * ($weighted_viremia[$sero_cnt][$other_region] / $annual_total_population[ $year ] ) * $seasonality;
				}

			}
	
			if ( $risk[$sero_cnt][$hmon][$region_pnt] == 0 ){
				# extincted virus revived.
				$risk[$sero_cnt][$hmon][$region_pnt] = $MINIMUM_RISK;
				$extinct_cnt[$sero_cnt][$year]++;
				print "e";
			}
		}

		my $sum_of_risk = 0;
		for ($region_pnt = 0; $region_pnt < $NUM_OF_REGIONS; $region_pnt++ ){
			$sum_of_risk += $risk[$sero_cnt][$hmon][$region_pnt];
		}
		my $average_risk = $sum_of_risk/$NUM_OF_REGIONS;
		
		$monthly_risk[$sero_cnt][$month] += $average_risk;

	}
}
######################################################
#  END OF Transmit_virus_to_environment              #
######################################################

######################################################
#  NaturalBirthDeath                                 #
######################################################
sub NaturalBirthDeath {

	my $human;
	my $sero_counter;

	if  ($NATURAL_AGE_STRUCTURE_FLAG > 0 ){

		&Natural_death;
			
		&Natural_birth;

	}
	else {

		for ($human=0; $human < $POPULATION; $human++){

			if ( $indiv_age_by_half_month[$human] >= $MAXIMUM_AGE_IN_YEAR*24 ){
				# a death is replaced by a birth
				
				$indiv_alive_or_dead[$human] = $ALIVE;
				
				for ($sero_counter = 1; $sero_counter <= $SEROTYPE; $sero_counter++ ){

					# epitope 1: protective epitope
					$indiv_protective_sero_specific_immunity_flag[ $sero_counter ][ $human ] = 0;
					$indiv_protective_epitope_sero_specific_inoculation_time_step[ $sero_counter ][ $human ] = 24*($SIM_YEAR+1000); # impossibly large number

				}

				$indiv_enhancement_by_wt_flag[ $human ] = 0;

				$indiv_enhancement_by_vaccine_flag[ $human ] = 0;

				$indiv_age_by_half_month[$human]=0;			#newborn

				$indiv_vaccinated_count[ $human ] = 0;
				$indiv_vaccinated_timestep[ $human ] = "";
			}
		}
	}
}
######################################################
#  END OF NaturalBirthDeath                          #
######################################################

######################################################
# Natural_death                                      #
######################################################
sub Natural_death {
	my $age_by_year;
	my $probability_to_die;
	my $age_class_cnt;

	for (my $human_cnt = 0; $human_cnt < $current_population_pointer; $human_cnt++ ){

		if ( $indiv_alive_or_dead[ $human_cnt ] == $ALIVE ){
			$age_by_year = int( $indiv_age_by_half_month[ $human_cnt ]/ 24);

			if ( $age_by_year > $MAXIMUM_AGE_IN_YEAR ){

				$indiv_alive_or_dead[ $human_cnt ] = $DEAD;

			}
			else {
				$age_class_cnt = int( $age_by_year / 5 ) + 1;
				if ( $age_class_cnt > $max_age_class ){
					$age_class_cnt = $max_age_class;
				}

				$probability_to_die = $age_specfic_annual_mortality_rate[ $age_class_cnt ] / 24;

				if ( rand() < $probability_to_die ){
					$indiv_alive_or_dead[ $human_cnt ] = $DEAD;
				}
			}
		}
	}
}
#####################################################################
# END OF Natural_death                                              #
#####################################################################

#####################################################################
# Natural_birth                                                     #
#####################################################################
sub Natural_birth {

	my $age_by_year;
	my $gender;
	my $age_class_cnt;
	my $dead_or_alive;
	my $number_of_child;
	my $sero_counter;
	my $probability_of_birth;
	
	for (my $human_cnt = 0; $human_cnt < $current_population_pointer; $human_cnt++ ){

		$age_by_year = int( $indiv_age_by_half_month[ $human_cnt ]/24) ;
		$dead_or_alive = $indiv_alive_or_dead[ $human_cnt ];
		$gender = $indiv_gender[ $human_cnt ];

		if ( $gender == $FEMALE ){
			if ( $dead_or_alive == $ALIVE ) {

				$age_class_cnt = 1 + int( $indiv_age_by_half_month[ $human_cnt ]/24/5 );
				
				if ( $age_class_cnt >= $min_fertile_age_class && $age_class_cnt <= $max_fertile_age_class ){

					$probability_of_birth = $age_specific_birthrate[ $age_class_cnt ]/24;

					if ( rand() < $probability_of_birth ){
					
						#birth

						for ($sero_counter = 1; $sero_counter <= $SEROTYPE; $sero_counter++ ){

							# epitope 1: protective epitope
							$indiv_protective_sero_specific_immunity_flag[ $sero_counter ][ $current_population_pointer ] = 0;
							$indiv_protective_epitope_sero_specific_inoculation_time_step[ $sero_counter ][ $current_population_pointer ] = 24*($SIM_YEAR+1000); # impossibly large number
						}

						$indiv_enhancement_by_wt_flag[ $current_population_pointer ] = 0;

						$indiv_enhancement_by_vaccine_flag[ $current_population_pointer ] = 0;

						$indiv_age_by_half_month[ $current_population_pointer ]=0;			#newborn

						$indiv_vaccinated_count[ $current_population_pointer ] = 0;
						$indiv_vaccinated_timestep[ $current_population_pointer ] = "";

						$indiv_alive_or_dead[ $current_population_pointer ] = $ALIVE;

						$indiv_region[ $current_population_pointer ] = $indiv_region[ $human_cnt ]; # inherit mother's region

						$indiv_gender[ $current_population_pointer ] = int(rand(2));

						$current_population_pointer++;
					}
				}
			}
		}
	}
}
#####################################################################
# END OF Natural_birth                                              #
#####################################################################


######################################################
# Transmit_virus_to_humans                           #
######################################################
sub Transmit_virus_to_humans {
	my $hmonth = $_[0];
	my $month = int ($hmonth/2);
	my $year = int ($hmonth/24);
	my $sero_cnt;
	my $human_cnt;
	my $age_year;
	my $int_age;
	my $age_class;
	my $experienced_protective_epitope_types;
	my $cross_protected_flag;
	my $sero_cnt_temp;
	my $elapsed_time;
	my $clinical_prob_prim;
	my $clinical_probability;
	my $my_region;
	my $temp_protective_period;
	my $temp_enhanced_viremia;
	my $temp_necessary_serotypes;

	for ($sero_cnt=1; $sero_cnt<= $SEROTYPE; $sero_cnt++ ){

		# initialize viremia
		
		for (my $region_pnt = 0; $region_pnt < $NUM_OF_REGIONS; $region_pnt++ ){
			$weighted_viremia[$sero_cnt][$region_pnt] = 0;	# number of viremic individuals weighted by ENHANCE factor.
		}
		
		for ($human_cnt = 0; $human_cnt < $current_population_pointer; $human_cnt++ ){

			if ( $indiv_alive_or_dead[ $human_cnt ] == $ALIVE ){
				$age_year = $indiv_age_by_half_month[$human_cnt] / 24;
				$int_age = int ($age_year);

				if ( $int_age < $CUTOFF_1 ){
					$age_class = 1;
				}
				elsif ( $int_age < $CUTOFF_2 ){
					$age_class = 2;
				}
				else {
					$age_class = 3;
				}

				$my_region = $indiv_region[ $human_cnt];
			
				if ( $risk[$sero_cnt][$hmonth][$my_region] >= rand() ){
					# biten by an infected mosquito

					$total_inoculation_number[$month]++;

					if ( $indiv_protective_sero_specific_immunity_flag[ $sero_cnt ][ $human_cnt] == 0 ){
						# no specific protective_epitope immunity.

						$experienced_protective_epitope_types = &Experienced_protective_epitope_types( $human_cnt );
						
						if ( $PDF_flag > 0 ){
							$temp_necessary_serotypes = int( rand(3) ) + 2; # 2 or 3 or 4
						}
						else {
							$temp_necessary_serotypes = $NECESSARY_SEROTYPES; # 2 or 3 or 4
						}

						if ( $experienced_protective_epitope_types >= $temp_necessary_serotypes ){

							#####################
							# totally immune    #
							#####################
							# acquire immunity to all serotypes
							for ($sero_cnt_temp = 1; $sero_cnt_temp <= $SEROTYPE; $sero_cnt_temp++ ){
								$indiv_protective_sero_specific_immunity_flag[ $sero_cnt_temp ][ $human_cnt] = 1;
								# immune to all serotypes
							}
							
						}
						else {
							$cross_protected_flag = 0;

							for ( $sero_cnt_temp = 1; $sero_cnt_temp <= $SEROTYPE; $sero_cnt_temp++ ){

								if ( $indiv_protective_sero_specific_immunity_flag[ $sero_cnt_temp ][ $human_cnt] > 0 ){
	
									$elapsed_time = $hmonth - $indiv_protective_epitope_sero_specific_inoculation_time_step[ $sero_cnt_temp ][ $human_cnt]; 
									# time steps which elapsed from the most recent inoculation or vaccination.

									if ( $PDF_flag > 0 ){
										$temp_protective_period = &randn( $PROTECTIVE_PERIOD, $PROTECTIVE_PERIOD/4 );
									}
									else {
										$temp_protective_period = $PROTECTIVE_PERIOD;
									}

									if ( $elapsed_time >= 0 && $elapsed_time <= $temp_protective_period ){
										# Cross-protection is effective.
										$cross_protected_flag = 1;
									}
								}
							}

							if ( $cross_protected_flag == 1 ){


								###################
								# CROSS_PROTECTED #
								###################
								# do nothing.
							}
							else {

								$number_of_viremic_individuals[$sero_cnt][$month]++;

								$weighted_viremia[$sero_cnt][$my_region]++;
								
								$num_of_viremic_individuals_in_timestep[$hmonth]++;


								if ( rand() < &Calc_DF_prob($age_year) ){
									$monthly_DF_number[ $month ]++;

									if ( $indiv_vaccinated_count[ $human_cnt ] > 0){
										$monthly_DF_number_in_vaccinated[ $month ]++;
									}
									else {
										$monthly_DF_number_in_non_vaccinated[ $month ]++;
									}

								}

								if ( $experienced_protective_epitope_types == 0 ){
					
									###################
									# NAIVE           #
									###################

									$clinical_prob_prim = $PRIMARY_DHF_PROB;

									if ( $clinical_prob_prim >rand() ){

										# DHF in primary infection manifested.

										print "p";

										if ( $PDF_flag > 0 ){
											$temp_enhanced_viremia = randn( $ENHANCED_VIREMIA, $ENHANCED_VIREMIA/4 );
										}
										else {
											$temp_enhanced_viremia = $ENHANCED_VIREMIA;
										}
										if ( $temp_enhanced_viremia < 1 ){
											$temp_enhanced_viremia = 1;
										}
										
										$weighted_viremia[$sero_cnt][$my_region] += ( $temp_enhanced_viremia - 1);
										# During DHF, the patient exerts enhanced transmissibility, defined by $ENHANCED_VIREMIA


										$annual_sum_of_age_of_DHF[ $year ] += $age_year;
										$DHF_number[$month][$sero_cnt]++;
										$total_DHF_number[ $month ]++;
										$primary_DHF_number[ $month ]++;

										$age_specific_DHF_number[$year][$age_class]++;

										if ( $indiv_vaccinated_count[ $human_cnt ] > 0){
											$DHF_number_in_vaccinated[ $month ]++;
											$DHF_number_in_vaccinated_serotype[$month][$sero_cnt]++;
										}
										else {
											$DHF_number_in_non_vaccinated[ $month ]++;
											$DHF_number_in_non_vaccinated_serotype[$month][$sero_cnt]++;
										}
									}
								}
								else {
									if ( $indiv_enhancement_by_wt_flag[ $human_cnt ] == 1 || $indiv_enhancement_by_vaccine_flag[ $human_cnt ] == 1){
										###################
										# DHF-predisposed #
										###################

										if ( $AGE_DEPEND == 0 ){ # no age dependency
											$clinical_probability = $SECONDARY_DHF_PROB;
										}
										elsif ( $AGE_DEPEND == 1 ){ # Severity decreases in parallel with age.
											$clinical_probability = ( (80 - $age_year )*0.06 + $age_year*0.02 )/80;
										}
										elsif ( $AGE_DEPEND == 2 ){ # Severity increases in parallel with age.
											$clinical_probability = ( ( 80 - $age_year )*0.02 + $age_year * 0.08 )/80;
										}
										if ( $clinical_probability >= 1 || $clinical_probability <= 0 ){
											print "secondary clinical probability wrong\n";
											exit;
										}

										if ( $clinical_probability >= rand()  ){

											print "h";

											if ( $PDF_flag > 0 ){
												$temp_enhanced_viremia = randn( $ENHANCED_VIREMIA, $ENHANCED_VIREMIA/4 );
											}
											else {
												$temp_enhanced_viremia = $ENHANCED_VIREMIA;
											}
											if ( $temp_enhanced_viremia < 1 ){
												$temp_enhanced_viremia = 1;
											}
											
											$weighted_viremia[$sero_cnt][$my_region] += ( $temp_enhanced_viremia - 1);
											$annual_sum_of_age_of_DHF[ $year ] += $age_year;

											$DHF_number[$month][$sero_cnt]++;
											$total_DHF_number[ $month ]++;

											$age_specific_DHF_number[$year][$age_class]++;


											if ( $indiv_vaccinated_count[ $human_cnt ] > 0 ){
												$DHF_number_in_vaccinated[ $month ]++;
												$DHF_number_in_vaccinated_serotype[$month][$sero_cnt]++;

												if ( $indiv_enhancement_by_vaccine_flag[ $human_cnt ] > 0 ){
													$annual_DHF_number_due_to_vaccine[ $year ]++;

													unless ( $indiv_enhancement_by_wt_flag[ $human_cnt ] > 0 ){
														$annual_DHF_number_due_to_vaccine_only[ $year ]++;
													}

												}

											}
											else {
												$DHF_number_in_non_vaccinated[ $month ]++;
												$DHF_number_in_non_vaccinated_serotype[$month][$sero_cnt]++;
											}
										
										}
									}
								}
							}
						}
					}

					$indiv_protective_sero_specific_immunity_flag[ $sero_cnt ][ $human_cnt] = 1;
					$indiv_protective_epitope_sero_specific_inoculation_time_step[ $sero_cnt ][ $human_cnt] = $hmonth;

					$indiv_enhancement_by_wt_flag[ $human_cnt] = 1;
				}
			}
		}
	}
}
######################################################
# END OF Transmit_virus_to_humans                    #
######################################################

######################################################
# Experienced_protective_epitope_types               #
######################################################
sub Experienced_protective_epitope_types {
	my $person = $_[0];
	my $experienced_protective_epitope_serotypes = 0;
	my $sero;

	for ($sero = 1; $sero <= $SEROTYPE; $sero++ ){
		if ( $indiv_protective_sero_specific_immunity_flag[ $sero ][ $person ] > 0 ){
			$experienced_protective_epitope_serotypes++;
		}
	}
	return($experienced_protective_epitope_serotypes);
}
######################################################
# END OF Experienced_infection                       #
######################################################

######################################################
# Calc_DF_prob                                       #
######################################################
sub Calc_DF_prob {
	my $age_by_year = $_[0];

	my $prob = 1/(1+1/exp(-3.44 + 0.177 * $age_by_year) );

	return ($prob);
}
######################################################
# END OF Calc_DF_prob                                #
######################################################

######################################################
# Calc_Effective_R                                   #
######################################################
sub Calc_Effective_R {
	my $hmonth = $_[ 0 ];
	my $year = int($hmonth / 24);
	
	if ( $num_of_viremic_individuals_in_timestep[$hmonth - 1] > 0 ){
		my $effective_R_of_timestep = $num_of_viremic_individuals_in_timestep[$hmonth]/$num_of_viremic_individuals_in_timestep[$hmonth - 1];
		$annual_sum_of_effective_R[ $year ] += $effective_R_of_timestep;
		$annual_cnt_of_effective_R[ $year ]++;
	}
}
######################################################
# END OF Calc_Effective_R                            #
######################################################

######################################################
# randn                                              #
######################################################
sub randn {
 my ($m, $sigma) = @_;
 my ($r1, $r2) = (rand(), rand());
 while ($r1 == 0) { $r1 = rand(); }
 return ($sigma * sqrt(-2 * log($r1)) * sin(2 * $pi * $r2)) + $m;
}
######################################################
# END OF randn                                       #
######################################################