##
## General note on program comments:
##
## Comments with "##" are standard to all my programs.
## Comments with "#" are descriptive of the specific program.
## Comments with "#." are lines of program code that can be used by deleting the "#.".
##

##
## General note on program structure:
## 
## It is generally best to look at the program elements in this order:
##  1. INTRODUCTION
##  2. OUTPUT
##  3. DATA
##  4. INITIALIZATION
##  5. MAIN PROGRAM
##  6. FUNCTIONS
##

##
## INTRODUCTION
##

# This program simulates the growth of an aphid population using a stage-structured model.

##
## DATA
##
## Import any data from files.
## Give values to any fixed parameters.
##

# define the model size and the number of days

# define the initial population of stage "i"

# define the fecundity of stage "i" 

# define the probability of development from stage "j" to stage "i=j+1"

# define the probability of survival in stage "i" without development 

# define the probability of the transition from oldest nymph to reproducing adult

##
## INITIALIZATION
##
## Set up any running variables that need to have a starting value.
## Create any needed data structures.
##

# define the data structures for the model, populations, total population, proportions, and daily growth rates

# enter the initial populations into "V"

# enter the "f" values into "M"

# enter the "p" values into "M" 

# enter the "s" values into "M" 

# enter "a" into "M"

##
## FUNCTIONS
##
## Define any functions that the main program requires.
## Calculations that need to repeated at different points of a program should be
##   coded as functions.
##

# "eigen(M)" is a built-in function that determines the eigenvalues and eigenvectors of a matrix "M"
#    "eigen(M)$values[1] is the largest eigenvalue
#    "eigen(M)$vectors[,1]" is an eigenvector for the largest eigenvalue


##
## MAIN PROGRAM
##
## Perform the computations, using functions as needed.
##

## FIND THE LONG-TERM GROWTH RATE AND AGE DISTRIBUTION

# Compute eigenvalues and eigenvectors using "eigen"

# "lambda1" is the largest eigenvalue

# "proportions" is the normalized eigenvector for the largest eigenvalue

# "r" is the continuous growth rate, which we need for any differential equation model


## RUN THE SIMULATION

# Compute "V"

# Compute "N" 			# "N[j]" is the total population at time "j-1"

# Compute "U"			# "U[,j]" is the column of proportions at time "j-1"

# Compute "R" 			# "R[j]" is the ratio of "N" at time "j" to "N" at time "j-1"


##
## OUTPUT
##
## Create any needed graphs.
## Display key results.
##

# determine the y-axis limit

# plot the data for daily growth rate

# connect the points in the growth rate graph

# plot the long-term growth rate as an asymptote

# print the output quantities

M			# matrix for the model
lambda1		# dominant eigenvalue
proportions		# proportions of stages for stable distribution
r			# growth rate for continuous model