Atmospheric column density

This function calculates atmospheric column density for one species in the solar direction (or X-axis of an arbitrary coordinate system)

A numerical integration of atmospheric column density for a given species at a point of interest (such as a position vector). Textbook examples are using solar zenith angle (polar angle counted from direction to the Sun) and altitude, but in practice it is much easier to use the position vectors (ecliptic coordinates) of the spacecraft or balloon in question. Either method can be used.

Inputs:

ecl_eph = [Xecl Yecl Zecl] matrix of ecliptical position vectors (in km) for the points of interest

IF method = 1, coord = [SZA alt] - the alt is point of interest, not necessarily same as in the density profile

N = [altitude density] matrix with a density profile (in cm^-3)

R_p = planet radius (in km)

hm = atmospheric top boundary altitude

vecstep = optional spatial integration step (in km), default is 1 km

the "method" variable (optional) is a switch for column density calculation method

if method = 0 or skipped, the ecliptical coordinates are used (i.e., default)

if method = 1, the SZA and Altitude are used

Output:

Ncol = column density (in cm^-2)

written by Oleg Shebanits, 2016

feel free to modify and adapt to your needs