(*10/19 need to fix string concatenation in the label of the graphs*)

(*These are routines to be used with bbseq/e*)


(* FExtract: Extracts a solution from bbseq/e output
	FResult: result of bbseq/e
	FPower: the iteration number to extract result from.
		If bbseq/e was run with stepdata which required	
		n iterations, with reporting every m steps (where m|n)
		then m <= FResult <= n
*)

FExtract[FResult_,FPower_] := Module[{r1,r2},
	r1=FResult/. u^FPower->B;
	r2=D[r1,B];
	psb00[[1]]=r2[[1]];
	psb00[[2]]=r2[[2]];
	Return[r2[[3]]];
	]

(*FEvalSoln: Evaluates how good a solution is using fsb.
	psb00: ps values
	eps00: eps rules
	r2: a solution 
*)

FEvalSoln[r2_, psb00_, eps00_] := Module[{r},
	ps={bz->psb00[[1]], b3->psb00[[2]], eps00[[1]], 
		eps00[[2]], eps00[[3]]};	
	prep[];
	r = r2//fsb;
	Return[r];
	]


(* IN POGRESS
This routine will take a solution, vary bz and compute how the solution divergesfrom 0. Outputs a 2 dimensional array.
r1 is a solution,  psb00 and eps00 are the same structures sent to bbseq/e*)

bzChange[r2_,psb00_,eps00_,steps_]:=Module[{bz0,bz1,stepsize, r, t1,ct, psb01},
	stepsize=1/steps;
	psb01=psb00;

	bz1=psb01[[1]];
	bz0=bz1-(IntegerPart[bz1*steps]*stepsize);
	
	r=Table[Null,{i,steps-1},{j,2}];
	ct=1;
	Label[o1];
	bz0=bz0+stepsize;
	psb01[[1]]=bz0;
	t1=FEvalSoln[r2,psb01,eps00];
	r[[ct]]={bz0,t1};
	ct=ct+1;
	If[bz0<=1-stepsize, Goto[o1]];

	Return[r];
	]
(* IN POGRESS
This routine will take a solution, vary b3 and compute how the solution diverges. Outputs a 2 dimensional array.
r1 is a solution,  psb00 and eps00 are the same structures sent to bbseq/e
*)

b3Change[r2_,psb00_,eps00_,steps_]:=Module[{b30,b31,stepsize, r, t1,ct, psb01},
	stepsize=1/steps;
	psb01=psb00;

	b31=psb01[[2]];
	b30=b31-(IntegerPart[b31*steps]*stepsize);
	If[b30<0,b30=0];
	r=Table[Null,{i,steps-1},{j,2}];
	ct=1;
	Label[o1];
	b30=b30+stepsize;
	psb01[[2]]=b30;
	t1=FEvalSoln[r2,psb01,eps00];
	r[[ct]]={b30,t1};
	ct=ct+1;
	If[b30<=1-stepsize, Goto[o1]];

	Return[r];
	]

(*Assumes ps is defined, parses ps and fills psb00, eps00 in the format required by bbseq/e *)
parsePs[psb01_]:=Module[{},
	psb00={0,0};
	eps00={0,0,0};
	psb00[[1]] = bz /.psb01[[1]];
	psb00[[2]] = b3 /.psb01[[2]];
	eps00[[1]]=psb01[[3]];
	eps00[[2]]=psb01[[4]];
	eps00[[3]]=psb01[[5]];	
	]

(*Run both b3Change and bzChange and generate graphs for each*)

bChange[r1_, psb00_, eps00_, steps_]:=Module[{sb3,sbz, lpba, lpb3},
	sb3=b3Change[r1,psb00,eps00,steps];
	lpb3=ListPlot[sb3, PlotLabel->"Variation with b3, bz" <> psb00[[1]], 
			AxesLabel->{"b3","fsb"}];
	sbz=bzChange[r1,psb00,eps00,steps];
	lpbz=ListPlot[sbz, PlotLabel->"Variation with bz, b3" <>  psb00[[2]], 
			AxesLabel->{"bz","fsb"}];	
	]

(*takes bbseq output, evaluates all the steps, and prints out fsb results
 so that you can choose the best one *)

FEvalAll[r2_,psb00_,eps00_,std00_]:=Module[{ct1, inc, st, tsol, fsbval},
	(*Step size*)
	ct1=0;
	inc=std00[[4]];
	(*endpoint*)
	st=std00[[3]];
	Label[o1];
	ct1=ct1+inc;
	tsol=FExtract[r2,ct1];
	fsbval= FEvalSoln[tsol,psb00,eps00];

	Print[ct1,": ", fsbval];
	If[ct1<st, Goto[o1]];

	Return[];
]

end=lastline