(* 4/11/95  modified to compute average value of absolute
	error between flux and accumulation curves *)
(* use this program to plot elem massbal from 1D fleet *)
(* modified so it will plot elem momementum balance also *)
(* ne is number of elements, elemlist contains element numbers 
you want to plot first element is element 1 *)

(* mass balance program *)
mass2[file1_String,ne_,elemlist_]:=Block[{i,graphlist={},dat,time,accum,netflux},
(*	SetOptions[ListPlot,DisplayFunction->Identity,Background->g10,Frame->True]; *)
	dat=ReadList[file1,Table[Number,{2 ne+2}]];
	time=take2[dat];time=time/3600.;
For[i=1,i<=Length[elemlist],i++,
	ii=elemlist[[i]];
	accum=taken[dat,2 ii+1]; accum=accum/10^6;
	netflux=taken[dat,2 ii+2]; netflux=netflux/10^6;
	error=Apply[Plus,Abs[accum-netflux]]/Length[accum];
	Print["Average of |mass_error| for element ", ii," = ",error];
	Print[];
(* comment out this Write if don't want to go to output file.
	NOTE!  Must open the output file (and close it) from
	the Mathematica session and must call it outfl *)
(*	Write[outfl,error]; *)
	gr1=ListPlot[xyjoin[time,netflux],PlotStyle->{d1,t1},FrameLabel->{"Time (hours)","Cumulative Mass (10^6)"},PlotLabel->StringJoin["Element #",ToString[ii]," (accum is solid)"]];
	gr2=ListPlot[xyjoin[time,accum],PlotStyle->{t1}];
	gr5=Show[gr1,gr2,Axes->False];
	disp[gr5];
(*	gr5=Show[gr5,Background->g10]; *)
	AppendTo[graphlist,gr5]
];
(* SetOptions[ListPlot,DisplayFunction:>$DisplayFunction]; *)
Return[graphlist]
	]

(*
(* momentum balance plots *)
mom2[ne_,elemlist_]:=Block[{i,graphlist={},dat,time,accum,netflux},
	SetOptions[ListPlot,DisplayFunction->Identity,Background->g10,Frame->True];
	dat=ReadList["~/newton/oneD/out6",Table[Number,{2 ne+2}]];
	time=take2[dat];time=time/3600.;
For[i=1,i<=Length[elemlist],i++,
	ii=elemlist[[i]];
	accum=taken[dat,2 ii+1]; accum=accum/10^6;
	netflux=taken[dat,2 ii+2]; netflux=netflux/10^6;
	gr1=ListPlot[xyjoin[time,netflux],PlotStyle->{d1,t1},FrameLabel->{"Time (hours)","Cumulative Momentum (10^6)"},PlotLabel->StringJoin["Element #",ToString[ii]," (accum is solid)"]];
	gr2=ListPlot[xyjoin[time,accum],PlotStyle->{t1}];
	gr5=Show[gr1,gr2,Axes->False];
	disp[gr5];
(*	gr5=Show[gr5,Background->g10]; *)
	AppendTo[graphlist,gr5]
];
SetOptions[ListPlot,DisplayFunction:>$DisplayFunction];
Return[graphlist]
	] *)