Problem II.1.1: Useful hints for numerical integration
==================================================================

The solar model you have now provides values for various
quantities of interest (density, temperature, etc.) at a 
finite set of N depth points. Using the radius as the independent
variable (second column in the model file), let 

r_j   j=,1...,N

represent the radius of the jth level, with similar notation
applying to other model quantities. 

In general, numerical integration seeks to replace the integral
operator (defined in terms of continuous variables) by a summation
(defined in terms of discrete quantities:

 R              _N
/               \
| f(r) dr  -->  / w_j f_j
/               -
 0               j=1

where w_i is some quadrature weight, f_j=f(r_j), and we
assumed that r_0=0 and r_N=R.

The simplest form of numerical integration is the TRAPEZOIDAL
RULE. This amounts to assuming a piecewise linear variation
of f(r) between succesive mesh points r_j, r_j+1.

 R              _N-1
/               \    1  
| f(r) dr  -->  /    -  (f_j+1 + f_j) * (r_j+1 - r_j) 
/               -    2  
 0               j=1

Note: the above formula is totally general and is applicable
to discretized variables for which the discretization grid
does not have constant spacing. For grids with constant spacing,
other formulae, involving higher order interpolation between
groups of mesh points, can be constructed.