/* Andrew D. Todd Simulation Model of the early English Population and Economy */ /* Population Simulation Section-- in progress. */ #DEFINE NUM_POPS 4 long huge big_arr[20640]; struct nutrit_total { float calories[NUM_POPS]; /*Kcal */ float protein[NUM_POPS]; /*grams */ float vitamins[NUM_POPS]; /* mg. Vit C */ }; /*-------------------------------------------------------------------*/ /* Interchange data area between population section and economy section. */ int time_slices = 10; long number_of_marriages_permitted[NUM_POPS]; struct nutrit_total Total_Normal_Food_Req; /*per year */ struct nutrit_total Total_Food_Avail_Harv; struct nutrit_total Total_Food_Avail_Now; struct nutrit_total Total_Food_Avail_New; /*-------------------------------------------------------------------*/ int age; /* 0 to 79 */ int last_age = 79; /*sexes*/ int sx_male = 1; int sx_female =2; int sx_first = 1 ; int sx_last = 2; /*marital status*/ int ms_single =1; int ms_married = 2; int ms_widowed = 3; int ms_first = 1; int ms_last = 3; /*sub-population*/ int sp_gentry = 1; int sp_rural = 2; int sp_urban = 3; int sp_first = 1; int sp_last = 3; /* program initialization variables */ long offset_mar_arr, t_arr_size; program_initialize() { offset_mar_arr = sp_last * ms_last * sx_last * (last_age+1); t_arr_size = offset_mar_arr + ((last_age +1) * (last_age +1) * sp_last); /* here is where a dynamic memory allocation would go */ } /*----------------------------------------------------------------------*/ /* */ /* The basic data arrays and their arithmetic functions */ /* */ /*----------------------------------------------------------------------*/ long Get_Pop( int age, int sex, int marital_stat, int sub_pop ) { long index_1_dim; index_1_dim = Array_Pop(age, sex, marital_stat, sub_pop); return(Get_Array_L(index_1_dim)); } /*----------------------------------------------------------------------*/ void Set_Pop( int age, int sex, int marital_stat, int sub_pop, long x ) { long index_1_dim; index_1_dim = Array_Pop(age, sex, marital_stat, sub_pop); Set_Array_L(index_1_dim, x); } /*----------------------------------------------------------------------*/ long Array_Pop( int age, int sex, int marital_stat, int sub_pop ) { long offset = 0; long size = 1; Array_Slice(&offset, &size, age, 0, last_age); Array_Slice(&offset, &size, sex, sx_first, sx_last); Array_Slice(&offset, &size, marital_stat, ms_first, ms_last); Array_Slice(&offset, &size, sub_pop, sp_first, sp_last); return(offset); } /*----------------------------------------------------------------------*/ long Add_Pop( int age, int sex, int marital_stat, int sub_pop, long x ) { long y; y=Get_Pop(age, sex, marital_stat, sub_pop); y+=x; Set_Pop(age, sex, marital_stat, sub_pop, y); return(y); } /* Returns New Value */ /*----------------------------------------------------------------------*/ long Red_by_Frac_Pop( int age, int sex, int marital_stat, int sub_pop, float fraction_to_reduce ) { long y, number_reduced; y=Get_Pop(age, sex, marital_stat, sub_pop); number_reduced=(long)((float)y)*fraction_to_reduce; y-=number_reduced; Set_Pop(age, sex, marital_stat, sub_pop, y); return(number_reduced); } /* returns the number the referenced cell is reduced by. */ /*----------------------------------------------------------------------*/ long Get_Marriages( int age_of_wife, int age_of_husb, int sub_pop ) { long index_1_dim; index_1_dim = Array_Marriages(age_of_wife, age_of_husb, sub_pop); return(Get_Array_L(index_1_dim)); } /*----------------------------------------------------------------------*/ void Set_Marriages( int age_of_wife, int age_of_husb, int sub_pop, long x ) { long index_1_dim; index_1_dim = Array_Marriages(age_of_wife, age_of_husb, sub_pop); Set_Array_L(index_1_dim, x); } /*----------------------------------------------------------------------*/ long Array_Marriages( int age_of_wife, int age_of_husb, int sub_pop ) { long offset = 0; long size = 1; Array_Slice(&offset, &size, age_of_wife, 0, last_age); Array_Slice(&offset, &size, age_of_husb, 0, last_age); Array_Slice(&offset, &size, sub_pop, sp_first, sp_last); offset += offset_mar_arr; return(offset); } /*----------------------------------------------------------------------*/ long Add_Marriages( int age_of_wife, int age_of_husb, int sub_pop, long x ) { long y; y=Get_Marriages(age_of_wife, age_of_husb, sub_pop); y+=x; Set_Marriages(age_of_wife, age_of_husb, sub_pop, y); return(y); } /*----------------------------------------------------------------------*/ long Red_by_Frac_Marriages( int age_of_wife, int age_of_husb, int sub_pop, float fraction_to_reduce ) { long y, number_reduced; y=Get_Marriages(age_of_wife, age_of_husb, sub_pop); number_reduced=(long)((float)y)*fraction_to_reduce; y-=number_reduced; Set_Marriages(age_of_wife, age_of_husb, sub_pop, y); return(number_reduced); } /* returns the number the referenced cell is reduced by. */ /*----------------------------------------------------------------------*/ long Array_Slice( long *accumulating_offset, long *accumulating_size, int index, int curr_low_bound, int curr_up_bound ) { int curr_dim_size; int curr_true_index; curr_dim_size = curr_up_bound - curr_low_bound + 1; curr_true_index = index - curr_low_bound; *accumulating_offset += ((long) curr_true_index) * (*accumulating_size); *accumulating_size *= (long) curr_dim_size; return(*accumulating_offset); } /*----------------------------------------------------------------------*/ long Get_Array_L(long index_1_dim) { return(big_arr[index_1_dim]); } /*----------------------------------------------------------------------*/ void Set_Array_L( long index_1_dim, long x ) { big_arr[index_1_dim] = x; } /*----------------------------------------------------------------------*/ /* */ /* Array Initializing Functions */ /* */ /*----------------------------------------------------------------------*/ void Advance_Ages() { int age; int age_of_wife; int age_of_husb; int sex; int marital_stat; int sub_pop; long x; /* the main population table */ for( sub_pop = sp_first; sub_pop <= sp_last; sub_pop++) for( marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(sex = sx_first ; sex <= sx_last; sex++) { for(age = last_age; age >= 1; age--) { x=Get_Pop(age-1, sex, marital_stat, sub_pop); Set_Pop(age, sex, marital_stat, sub_pop, x); } Set_Pop(0, sex, marital_stat, sub_pop, 0); } /* Initialize age 0 */ /* and a supplemmentary table expressing the relative ages of spouses */ for( sub_pop = sp_first; sub_pop <= sp_last; sub_pop++) for(age_of_wife = last_age; age_of_wife >= 0; age_of_wife--) for(age_of_husb = last_age; age_of_husb >= 0; age_of_husb--) { if((age_of_wife == 0) || (age_of_husb == 0)) x = 0; else x=Get_Marriages( age_of_wife-1, age_of_husb-1, sub_pop ); Set_Marriages( age_of_wife, age_of_husb, sub_pop, x ); } } /* copies the population figures for every age into the next older age, simulating the passage of a year, initializes the new entries by setting them to zero */ /*----------------------------------------------------------------------*/ void Initialize_Pop_Tables(int sub_pop, long target_pop) { int age; int age_of_wife; int age_of_husb; int sex; int marital_stat; long entry_target_pop; long number_of_entries; /* the main population table */ number_of_entries = (last_age+1) * sx_last; entry_target_pop = (long) ( ( (float) target_pop) / ( (float) number_of_entries) ); for(marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(sex = sx_first ; sex <= sx_last; sex++) for(age = last_age; age >= 0; age--) { if(marital_stat == ms_single) Set_Pop(age, sex, marital_stat, sub_pop, entry_target_pop); else Set_Pop(age, sex, marital_stat, sub_pop, 0); } /* and a supplemmentary table expressing the relative ages of spouses */ for(age_of_wife = last_age; age_of_wife >= 0; age_of_wife--) for(age_of_husb = last_age; age_of_husb >= 0; age_of_husb--) Set_Marriages(age_of_wife, age_of_husb, sub_pop, 0); } /* sets entries in both tables to reflect a specified population number for the specified sub-population, all single, with equal numbers of both sexes, evenly distributed for age. */ /*----------------------------------------------------------------------*/ void Normalize_Pop_Tables(int sub_pop, long target_pop) { int age; int age_of_wife; int age_of_husb; int sex; int marital_stat; long x,y; float fraction_to_reduce; /* the main population table */ y=Total_Pop(sub_pop); fraction_to_reduce=(float)(target_pop-y)/((float)y); for( marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(sex = sx_first ; sex <= sx_last; sex++) for(age = last_age; age >= 0; age--) Red_by_Frac_Pop(age, sex, marital_stat, sub_pop, fraction_to_reduce); /* and a supplementary table expressing the relative ages of spouses */ y=Total_Marriages(sub_pop); fraction_to_reduce=(float)(target_pop-y)/((float)y); for(age_of_wife = last_age; age_of_wife >= 0; age_of_wife--) for(age_of_husb = last_age; age_of_husb >= 0; age_of_husb--) Red_by_Frac_Marriages( age_of_wife, age_of_husb, sub_pop, fraction_to_reduce ); } /* reduces every entry in the population tables for a given sub-population proportionately to reduce the total population to a stated figure. Should work even if target population is greater than actual population. */ /*----------------------------------------------------------------------*/ long Total_Pop(int sub_pop) { int age; int sex; int marital_stat; long x; x=0; for( marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(sex = sx_first ; sex <= sx_last; sex++) for(age = last_age; age >= 0; age--) x+=Get_Pop(age, sex, marital_stat, sub_pop); return(x); } /*----------------------------------------------------------------------*/ long Total_Marriages_By_Age( int sex, int age_of_sex, int sub_pop ) { long total; int age; total = 0; if(sex == sx_male) for(age = 0; age <= last_age; age++) total += Get_Marriages(age, age_of_sex, sub_pop); else for(age = 0; age <= last_age; age++) total += Get_Marriages(age_of_sex, age, sub_pop); return(total); } /*----------------------------------------------------------------------*/ long Total_Marriages(int sub_pop) { long total; int f_age, m_age; total = 0; for(f_age = 0; f_age <= last_age; f_age++) for(m_age = 0; m_age <= last_age; m_age++) total += Get_Marriages(f_age, m_age, sub_pop); return(total); } /*----------------------------------------------------------------------*/ /* */ /* Compute Food Requirement */ /* */ /*----------------------------------------------------------------------*/ void Get_Normal_Food_Req() { int age, sex, marital_stat, sub_pop; for( sub_pop = sp_first; sub_pop <= sp_last; sub_pop++) { Total_Normal_Food_Req.calories[sub_pop] = 0.0; Total_Normal_Food_Req.protein[sub_pop] = 0.0; Total_Normal_Food_Req.vitamins[sub_pop] = 0.0; for( marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(sex = sx_male ; sex <= sx_female; sex++) for(age = 0; age <= last_age; age++) { Total_Normal_Food_Req.calories[sub_pop] += Indiv_Normal_Food_Req_cal(age,sex) * (float) Get_Pop(age, sex, marital_stat, sub_pop); Total_Normal_Food_Req.protein[sub_pop] += Indiv_Normal_Food_Req_pro(age,sex) * (float) Get_Pop(age, sex, marital_stat, sub_pop); Total_Normal_Food_Req.vitamins[sub_pop] += Indiv_Normal_Food_Req_vit(age,sex) * (float) Get_Pop(age, sex, marital_stat, sub_pop); } } } /* Computes the ideal food requirement (per year) for the various sub-populations, broken down into: calories (esp fr. grain and potatoes), protein (esp fr. meat) and vitamins. called by the model advance */ /*-------------------------------------------------------------------------*/ float Indiv_Normal_Food_Req_cal(int age, int sex) { float desired_wt_kg, daily_kcal_per_kg; if( age == 0 ) desired_wt_kg = 7.5; else if( age < 4 ) desired_wt_kg = 13.0; else if( age < 7 ) desired_wt_kg = 20.0; else if( age < 11 ) desired_wt_kg = 28.0; else if(sex == sx_male) { if(age < 15 ) desired_wt_kg = 45.0; else if( age < 19 ) desired_wt_kg = 66.0; else desired_wt_kg = 70.0; } else if(sex == sx_female) { if(age < 15 ) desired_wt_kg = 46.0; else desired_wt_kg = 70.0; } /* always */ if( age == 0 ) daily_kcal_per_kg = 110.0; else if( age < 4 ) daily_kcal_per_kg = 100.0; else if( age < 7 ) daily_kcal_per_kg = 90.0; else if( age < 10 ) daily_kcal_per_kg = 80.0; else if( age < 13 ) daily_kcal_per_kg = 70.0; else if( age < 16 ) daily_kcal_per_kg = 60.0; else if( age < 20 ) daily_kcal_per_kg = 50.0; else daily_kcal_per_kg = 40.0; /*always*/ return( daily_kcal_per_kg * desired_wt_kg * 365.25); } /* returns Kcal per year, Merck, pp. 876-77, 1847 */ /*-------------------------------------------------------------------------*/ float Indiv_Normal_Food_Req_pro(int age, int sex) { if(age == 0) return(38.25 * 365.25); /* based on body weight, assume average 7.5 kg., * 2.1 gm per kg. 15.75 for the infant, 22.5 for the pregnant mother (30 prorated to 9 months) */ else if(age < 4 ) return(23.0 * 365.25); else if(age < 7 ) return(30.0 * 365.25); else if(age < 11 ) return(34.0 * 365.25); else if(sex == sx_male) { if(age < 15 ) return( 45.0 * 365.25); else return( 56.0 * 365.25 ); } else if(sex == sx_female) { if(age < 19 ) return( 46.0 * 365.25); else return( 44.0 * 365.25 ); } } /* returns grams per year, Merck, p. 876-77 */ /*-------------------------------------------------------------------------*/ float Indiv_Normal_Food_Req_vit(int age, int sex) { if(age == 0) return(50.0 * 365.25); /* 35 for the infant, 15 for the pregnant mother (20 prorated to 9 months) */ else if(age < 11 ) return(45.0 * 365.25); else if(sex == sx_male) { if(age < 15 ) return( 50.0 * 365.25); else return( 60.0 * 365.25 ); } else if(sex == sx_female) { if(age < 15 ) return( 50.0 * 365.25); else return( 60.0 * 365.25 ); } } /*returns mg. Vitamin C per year, Merck, p. 876-77 */ /*----------------------------------------------------------------------*/ /* */ /* Food Distribution Computations */ /* */ /*----------------------------------------------------------------------*/ void Draw_Down_Harvest(int slice_num){} /* takes food from Total_Food_Avail_Harv, Total_Food_Avail_New and puts in in Total_Food_Avail_Now */ /*----------------------------------------------------------------------*/ /* */ /* Mortality Computations */ /* */ /*----------------------------------------------------------------------*/ float Mortality( int age, int sex, int marital_stat, float General_Mortality_Level ) { int lbound, hbound, i; float Interpolation_fraction, General_Life_Expectancy; General_Life_Expectancy = 1/General_Mortality_Level; /* find the entries in the life expectancy table which straddle this. */ lbound = 0; hbound = 9; for(i=0; i < 10; i++) { if (mort_table_e_sub_0_entry(i) < General_Life_Expectancy) lbound = i; else { hbound = i; break; } } /* and interpolate between them */ Interpolation_fraction = (General_Life_Expectancy - mort_table_e_sub_0_entry(lbound) ) / (mort_table_e_sub_0_entry(hbound) - mort_table_e_sub_0_entry(lbound) ); return( Mortality_col(age, lbound) + Interpolation_fraction * ( Mortality_col(age, hbound) - Mortality_col(age, lbound) ) ); } float Mortality_col(int age, int col) { float table_col_mort; if(age = 0) table_col_mort = mort_table_mort_entry(0, col); else if(age < 5) table_col_mort = mort_table_mort_entry(1, col)/4; else if(age < 90) table_col_mort = mort_table_mort_entry( ((age/5)+1), col )/5; /* AGE is an integer, and is rounded down to the nearest whole number */ else table_col_mort = mort_table_mort_entry(19, col)/5; return(table_col_mort/1000); } /* Mortality Table Levels: 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 Life Expectancy: */ float mort_table_e_sub_0_arr[10] = {23.54, 25.98, 28.42, 30.86, 33.31, 35.77, 38.23, 40.7, 43.18, 45.66}; float mort_table_mort_arr[200] = { /* 0 1 2 3 4 5 6 7 8 9 */ /*age*/ /* 0*/ 279.2, 255.8, 234.4, 214.7, 196.6, 179.8, 164.1, 149.5, 135.8, 123.0,/*0*/ /* 1*/ 265.7, 241.8, 220.1, 200.1, 181.6, 164.5, 148.6, 133.7, 119.7, 106.6,/*1*/ /* 5*/ 91.9, 83.7, 76.2, 69.4, 63.0, 57.2, 51.7, 46.6, 41.8, 37.3,/*2*/ /*10*/ 49.3, 45.0, 41.1, 37.5, 34.2, 31.1, 28.2, 25.5, 23.0, 20.7,/*3*/ /*15*/ 58.3, 53.5, 49.2, 45.2, 41.5, 38.0, 34.8, 31.9, 29.1, 26.5,/*4*/ /*20*/ 77.6, 71.4, 65.7, 60.6, 55.8, 51.3, 47.2, 43.3, 39.7, 36.3,/*5*/ /*25*/ 86.3, 79.3, 73.0, 67.2, 61.8, 56.8, 52.2, 47.8, 43.8, 39.9,/*6*/ /*30*/ 95.4, 87.7, 80.6, 74.1, 68.1, 62.6, 57.4, 52.6, 48.0, 43.8,/*7*/ /*35*/ 106.5, 97.8, 89.8, 82.6, 75.9, 69.6, 63.8, 58.4, 53.3, 48.6,/*8*/ /*40*/ 118.4, 109.0, 100.4, 92.5, 85.2, 78.4, 72.1, 66.2, 60.7, 55.5,/*9*/ /*45*/ 135.0, 124.5, 114.9, 106.1, 98.0, 90.4, 83.4, 76.8, 70.7, 64.9,/*0*/ /*50*/ 158.2, 146.7, 136.2, 126.5, 117.6, 109.3, 101.6, 94.3, 87.6, 81.2,/*1*/ /*55*/ 205.8, 191.0, 177.4, 164.9, 153.3, 142.6, 132.6, 123.3, 114.5, 106.0,/*2*/ /*60*/ 275.2, 256.4, 239.2, 223.4, 208.8, 195.2, 182.6, 170.7, 159.7, 149.0,/*3*/ /*65*/ 377.8, 353.4, 331.2, 310.7, 291.9, 274.4, 258.1, 242.8, 228.6, 214.9,/*4*/ /*70*/ 527.2, 495.4, 466.3, 439.6, 414.9, 392.0, 370.7, 350.7, 332.0, 313.8,/*5*/ /*75*/ 682.4, 646.1, 613.0, 582.5, 554.3, 528.1, 503.7, 480.9, 459.5, 438.4,/*6*/ /*80*/ 790.1, 759.3, 731.2, 705.3, 681.3, 659.1, 638.3, 618.9, 600.8, 582.8,/*7*/ /*85*/ 888.4, 867.1, 847.6, 829.6, 813.0, 797.6, 783.2, 769.8, 757.2, 744.8,/*8*/ /*90*/ 944.8, 932.9, 922.1, 912.1, 902.9, 894.3, 886.3, 878.8, 871.8, 864.9 /*9*/ }; /* After Wrigley and Schoefield, p. 714 These entries are the mortality from the age associated with the row to the age associated with the next row, in deaths/1000. */ float mort_table_e_sub_0_arr(int i) { return( (float) mort_table_e_sub_0_arr[i] ); } float mort_table_mort_entry(int row, int col) { return( mort_table_mort_arr[row*10+col]); } /* This table is the means whereby a general population mortality is translated into age mortality */ /*----------------------------------------------------------------------*/ float Population_Mortality_Level(int sub_pop) { float Basal_Mortality = 0.014; float Hunger_Coefficient = 4.0 ; /* a guess-- complete lost of food results in starvation in 3 months */ float Crowding_Mortality; if( sub_pop = sp_gentry) Crowding_Mortality = 0.0; else if( sub_pop = sp_rural) Crowding_Mortality = 0.1; else if( sub_pop = sp_urban) Crowding_Mortality = 0.05; /* guesses, an urban prole can expect to get a fatal communicable disease once in twenty years, a peasant, once in a hundred years */ /*always*/ return( Basal_Mortality + Hunger_Coefficient * (1.0 - Compute_Fract_Food_Required(sub_pop)) + Crowding_Mortality ); } /* This determines the general mortality level for a population. Environmental effects will be introduced here. Provisionally mortality will be that of modern men and hunter-gatherers (about 1/70) with deductions for food shortfall and population density. */ /*----------------------------------------------------------------------*/ float Compute_Fract_Food_Required(int sub_pop) { float Frac_Calories_Avail, Frac_Protein_Avail, Frac_Vitamins_Avail, Calorie_Deficit; float Calorie_BaseLine = 0.5; /* The highest value for calories alone */ float Protein_BaseLine = 0.8; /* ditto for calories and protein */ /* Wild Guesses*/ /* if sufficient calories are not available, metabolize protein, at the rate of 4 Kcal/gm (Merck, p. 876) */ Calorie_Deficit = ( Total_Normal_Food_Req.calories[sub_pop] / (float) time_slices ) - Total_Food_Avail_Now.calories[sub_pop]; if( Calorie_Deficit > 0.0 ) /* convert some protein */ { if( Total_Food_Avail_Now.protein[sub_pop] * 4.0 /* Kcal per gm*/ >= Calorie_Deficit ) { Total_Food_Avail_Now.calories[sub_pop] = Total_Normal_Food_Req.calories[sub_pop] / (float) time_slices; Total_Food_Avail_Now.protein[sub_pop] -= Calorie_Deficit / 4.0; Frac_Calories_Avail = 1.0; } else /* not sufficient protein */ { Total_Food_Avail_Now.calories[sub_pop] += Total_Food_Avail_Now.protein[sub_pop] * 4; /* Kcal per gm*/ Total_Food_Avail_Now.protein[sub_pop] = 0.0; Frac_Calories_Avail = Total_Food_Avail_Now.calories[sub_pop] * time_slices / Total_Normal_Food_Req.calories[sub_pop]; } } else { Frac_Calories_Avail = 1.0; } /* Always */ /*if there is still a calorie deficit*/ if( Frac_Calories_Avail < 1.0) { Total_Food_Avail_Now.vitamins[sub_pop] = 0.0; Frac_Vitamins_Avail = 0.0; /* The tiny quantities of vitamins are used up without any meaningful gain in calories */ } else /* supposing sufficient calories */ { Frac_Vitamins_Avail = Total_Food_Avail_Now.vitamins[sub_pop] * time_slices / Total_Normal_Food_Req.vitamins[sub_pop]; Frac_Protein_Avail = Total_Food_Avail_Now.protein[sub_pop] * time_slices / Total_Normal_Food_Req.protein[sub_pop]; } if(Frac_Calories_Avail < 1.0) return( Frac_Calories_Avail * Calorie_BaseLine ); else if(Frac_Protein_Avail < 1.0) return( Calorie_BaseLine + ( Frac_Protein_Avail * (Protein_BaseLine-Calorie_BaseLine) ) ); else /*if Frac_Protein_Avail >= 1.0*/ return (Protein_BaseLine +(Frac_Vitamins_Avail*(1.0 - Protein_BaseLine)) ); } /* returns a number between 0.00 (immediate starvation) and 1.00 (complete nutritional health) */ /*----------------------------------------------------------------------*/ void Apply_Mortality() { int age, sex, marital_stat, sub_pop; float mortality_this_time; float pop_mort_this_time; long deaths_this_time; for( sub_pop = sp_first; sub_pop <= sp_last; sub_pop++) { pop_mort_this_time = Population_Mortality_Level(sub_pop); for( marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(sex = sx_male ; sex <= sx_female; sex++) for(age = 0; age <= last_age; age++) { mortality_this_time=Mortality(age, sex, marital_stat, pop_mort_this_time ) / time_slices; deaths_this_time=Red_by_Frac_Pop(age, sex, marital_stat, sub_pop, mortality_this_time); if(marital_stat == ms_married) { if(sex ==sx_male) Widow_Wives(age, sub_pop, deaths_this_time); else Widow_Husbs(age, sub_pop, deaths_this_time); } } } } /* This obtains the mortality rate for every group and adjusts the population accordingly. It also invokes the routine to simulate widowhood. */ /*----------------------------------------------------------------------*/ Widow_Wives(int age_of_husb, int sub_pop, long y) { int age_of_wife; float Fraction_Mortality_Of_Husbs, Deaths_of_Husbs; Fraction_Mortality_Of_Husbs = F_Div_L( y, Total_Marriages_By_Age(sx_male, age_of_husb, sub_pop) ); /* reduce each entry of the right age by that percentage */ for(age_of_wife = 0; age_of_wife <= last_age; age_of_wife++) { Deaths_of_Husbs = Red_Frac_Marriages( age_of_wife, age_of_husb, sub_pop, Fraction_Mortality_Of_Husbs ); Widow_Pop( age_of_wife, sx_female, sub_pop, Deaths_of_Husbs ); } } /* removes the appropriate number of husbands from the appropriate diagonal of the husbands table, and shifts the wives from married to widowed in the main table */ /*----------------------------------------------------------------------*/ Widow_Husbs(int age_of_wife, int sub_pop, long y) { int age_of_husb; float Fraction_Mortality_Of_Wives, Deaths_of_Wives; Fraction_Mortality_Of_Wives = F_Div_L( y, Total_Marriages_By_Age(sx_female, age_of_wife, sub_pop) ); /* reduce each entry of the right age by that percentage */ for(age_of_husb = 0; age_of_husb <= last_age; age_of_husb++) { Deaths_of_Wives = Red_Frac_Marriages( age_of_wife, age_of_husb, sub_pop, Fraction_Mortality_Of_Wives ); Widow_Pop( age_of_husb, sx_male, sub_pop, Deaths_of_Wives ); } } /* removes the appropriate number of husbands from the appropriate line of the husbands table, and shifts them from married to widowed in the main table. */ /*----------------------------------------------------------------------*/ float F_Div_L(x,y) { float frac; frac= ( (float) x) / ( (float) y); return(frac); } /*----------------------------------------------------------------------*/ Widow_Pop( int age, int sex, int sub_pop, long N_Widowed ) { Add_Pop(age, sex, ms_widowed, sub_pop, N_Widowed); Add_Pop(age, sex, ms_married, sub_pop, (-N_Widowed)); } /*----------------------------------------------------------------------*/ Move_Pop( int src_age, int src_sex, int src_marital_stat, int src_sub_pop, int dst_age, int dst_sex, int dst_marital_stat, int dst_sub_pop, long N_moved ) { Add_Pop(src_age, src_sex, src_marital_stat, src_sub_pop, (-N_moved)); Add_Pop(dst_age, dst_sex, dst_marital_stat, dst_sub_pop, N_moved); } /*----------------------------------------------------------------------*/ /* */ /* Fertility Computation */ /* */ /*----------------------------------------------------------------------*/ float Fertility(int age_of_mother, float General_Population_Fertility_Level) { float age_fertility, gross_fertility, net_fertility; /*--------- look age up in fertility table -------------------------*/ if( age < 20 ) age_fertility = 0.380; /* assumed to be equal to value for age 20-24 */ else if( age < 25 ) age_fertility = 0.380; /* W&S, Table 7.25 */ else if( age < 30 ) age_fertility = 0.341; /* p. 254 */ else if( age < 35 ) age_fertility = 0.288; /* giving fertillity */ else if( age < 40 ) age_fertility = 0.236; /* for ages 20-44 */ else if( age < 45 ) age_fertility = 0.121; /* */ else age_fertility = 0.0; /* assumed to be 0 */ /*-------- then correct for overall fertillity level ---------------*/ gross_fertility = age_fertility * General_Population_Fertility_Level; if( gross_fertility >= 1 ) net_fertility = 1; else net_fertility = gross_fertility; return(net_fertility); } /* This function takes W&S's mean values for various ages as a baseline, and multiplies them by the General Population Fertility Level. The result is capped at one child per woman per year. */ /*----------------------------------------------------------------------*/ float Population_Fertility_Level(int sub_pop) { return(1); } /* This function will reflect environmental effects, which is to say nutrition. However, there seems little case for environmental effects at this time, so set it to 1. */ /*----------------------------------------------------------------------*/ Apply_Fertility() { int age; int sex; int marital_stat; int sub_pop; long total, subtotal, N_Potential_Mothers; float Effective_Fertility; for( sub_pop = sp_first; sub_pop <= sp_last; sub_pop++) { total=0; for( marital_stat = ms_first; marital_stat <= ms_last; marital_stat++) for(age = 0; age <= last_age; age++) { N_Potential_Mothers = Get_Pop(age, sx_female, marital_stat, sub_pop); Effective_Fertility = Fertility(age, Population_Fertility_Level(sub_pop)) * Fem_Sex_Particip_Rate(age, marital_stat, sub_pop) / time_slices; subtotal = (long)( (float) N_Potential_Mothers * Effective_Fertility); total += subtotal; } Set_Pop(0, sx_male, ms_single, sub_pop, (total/2)); Set_Pop(0, sx_female, ms_single, sub_pop, (total/2)); } } /* This obtains the birth rates for every population, takes account of the number of married women, plus an allowance for the unmarried ones, computes the birth totals and inserts them into the population table. */ /*----------------------------------------------------------------------*/ float Fem_Sex_Particip_Rate(int age, int marital_stat, int sub_pop) { if(marital_stat == ms_married) return(1.0); else return(0.0); } /* Here is where we would insert the bastardy rate */ /*----------------------------------------------------------------------*/ /* */ /* Marriage Computation */ /* */ /*----------------------------------------------------------------------*/ Apply_Marriage() { int sub_pop; for(sub_pop = sp_first; sub_pop <= sp_last; sub_pop++) Marry_Off_Sub_Pop(sub_pop, number_of_marriages_permitted[sub_pop]); } /* */ /*-------------------------------------------------------------------------*/ long Marry_Off_Sub_Pop( int sub_pop, long total_number_permitted ) { long total_in_existence, number_performed, number_to_perform, number_still_to_perform, total_number_performed; if( total_number_permitted > (total_in_existence = Total_Marriages(sub_pop)) ) { number_to_perform = total_number_permitted - total_in_existence; number_still_to_perform = number_to_perform; total_number_performed = 0; number_performed = 1; while( (number_performed > 0) && (number_still_to_perform > 0) ) { number_performed = Perform_Some_Marriages(sub_pop, number_still_to_perform); number_still_to_perform -= number_performed; total_number_performed += number_performed; } return(total_number_performed); } else return(0); } /* Counts the number of marriages in being, and performs additional marriages to bring the number up to the total permitted. Returns the number actually performed. */ /*-------------------------------------------------------------------------*/ long Perform_Some_Marriages( int sub_pop, long additional_number_permitted ) { int f_age, f_marital_stat, m_age, m_marital_stat; long number_found, number_performed; Find_Best_Marital_Pair( &f_age, &f_marital_stat, &m_age, &m_marital_stat, sub_pop, &number_found ); if(number_found <= additional_number_permitted) number_performed = number_found; else number_performed = additional_number_permitted; Move_Pop( f_age, sx_female, f_marital_stat, sub_pop, f_age, sx_female, ms_married, sub_pop, number_performed ); Move_Pop( m_age, sx_male, m_marital_stat, sub_pop, m_age, sx_male, ms_married, sub_pop, number_performed ); Add_Marriages(f_age, m_age, sub_pop, number_performed); return(number_performed); } /* Returns number actually performed. If number_permited is greater than zero, return of zero indicates that there is no one left to marry. */ /*-----------------------------------------------------------------------*/ int Find_Best_Marital_Pair( int *f_age, int *f_marital_stat, int *m_age, int *m_marital_stat, int sub_pop, long *number_found ) { long f_number_found, m_number_found; Find_Best_Marital_Party( f_age, sx_female, f_marital_stat, sub_pop, &f_number_found ); Find_Best_Marital_Party( m_age, sx_male, m_marital_stat, sub_pop, &m_number_found ); if( ( f_number_found <= 0 ) || ( m_number_found <= 0) ) { *number_found = 0; return(0); } else if( f_number_found <= m_number_found ) { *number_found = f_number_found; return(1); } else { *number_found = m_number_found; return(1); } } /*----------------------------------------------------------------------*/ int Find_Best_Marital_Party( int *age, int sex, int *marital_stat, int sub_pop, long *number_found ) { int t_age, t_marital_stat; float Peak_Marital_Desirability, Curr_Marital_Desirability; Peak_Marital_Desirability = 0; for( t_marital_stat = ms_first; t_marital_stat <= ms_last; t_marital_stat++) for(t_age = 0; t_age <= last_age; t_age++) { if(t_marital_stat != ms_married) { Curr_Marital_Desirability = Marital_Desirability(t_age, sex, t_marital_stat, sub_pop); if( ( Curr_Marital_Desirability > Peak_Marital_Desirability ) && ( Get_Pop(t_age, sex, t_marital_stat, sub_pop) > 0 ) ) { *age=t_age; *marital_stat=t_marital_stat; Peak_Marital_Desirability = Curr_Marital_Desirability; } } } if( Peak_Marital_Desirability > 0 ) { *number_found = Get_Pop(*age, sex, *marital_stat, sub_pop); return(1); } else { number_found = 0; return(0); } } /*----------------------------------------------------------------------*/ float Marital_Desirability(age,sex,marital_stat,sub_pop) {} /* I am contemplating using a Weibull function: f(x) = alpha *(beta**(-alpha)) *(x**(alpha-1)) *(E**(-( (x/beta) ** alpha) )) if x > 0, else f(x) = 0; As you can see, it's a doozy, but Law and Kelton have some plots of it, which are about the right shape I want, with alpha of about 2-3. Ref: Averil M. Law and W. David Kelton, Simulation Modeling and Analysis, McGraw-Hill, New York, 1982, pp. 161-63 */ /*----------------------------------------------------------------------*/ /* */ /* Population Model Mechanisms */ /* */ /*----------------------------------------------------------------------*/ Advance_Population_Model() { int i; Advance_Ages(); for(i = 1; i >= time_slices; i++) { Get_Normal_Food_Req(); Draw_Down_Harvest(i); Apply_Fertility(); Apply_Mortality(); /* (and widowhood) */ Apply_Marriage(); } Get_Normal_Food_Req(); } /*----------------------------------------------------------------------*/ /* I propose to use an "adam and eve" start */ main(){} simulation_pass(){} /* sequence is 1. advance simulated time 2. advance popuation model 3. simulate the economy 4. and do it again */