function moment,data;

; procedure calculates the zero, first, second,and third moment for

; a breakthrough curve. The data for the breakthrough curve is in data.

; data consists of 2 columns a date_time column and a conc column

;info, data

;dt_print, data.date_time(0)

date_time=data.date_time

date_time=date_time.julian

date_time=date_time-date_time(0)

conc=data.conc

;perform two integrations first the area under the breakthrough

; curve and then the element times the moment arm */

last_row = N_Elements(data.(1))

sum_0 = 0.; /* zero moment */

for i = 0, last_row - 2 do begin;

dt = (date_time(i+1)) - (date_time(i));

c_bar = ((conc(i+1)) + (conc(i)))/2;

sum_0 = sum_0 + (dt * c_bar);

endfor;

 

sum_a = 0.;

sum_b = 0.;

sum_c = 0.;

for i = 0, last_row -2 do begin;

t1 = date_time(i) ;

t2 = date_time(i+1) ;

c1 = conc(i) ;

c2 = conc(i+1) ;

m_1 = float((((c1*t1)+(c2*t2))*(t2-t1) + (c1+c2)*((t2*t2)-(t1*t1)))/6.)

m_2 = float((((c1*t1*t1+c2*t2*t2)*(t2-t1)+(c1*t1+c2*t2)*(t2*t2-t1*t1)$

+(c1+c2)*(t2*t2*t2-t1*t1*t1)))/12.)

m_3 =float((((c1*t1^3+c2*t2^3)*(t2-t1)+(c1*t1*t1+c2*t2*t2)*(t2*t2-t1*t1)$

+(c1*t1+c2*t2)*(t2^3-t1^3)+(c1-c2)*(t2^4-t1^4)))/20.)

sum_a = sum_a + m_1

sum_b = sum_b + m_2

sum_c = sum_c + m_3

endfor;

mom=fltarr(4)

mom(0)=sum_0

mom(1)=sum_a

mom(2)=sum_b

mom(3)=sum_c

;print, mom

finish:

return, mom

end