Print["MBnum v.0.1, last modified: 18.06.09"];


MBnum::usage = "MBnum[repr_, epsilon_, points_, powers_, loop_]"


Options[MBnum]={
	RulesConstrain->60,
	ShowMBrep->False,
	ShowPositions->True,
	ShowRules->False,
	ShowComments->True,
	ShowMBcontinue->False,
	ShowMBexpand->False,
	ShowMBintegrate->False,
	ShowETArules->False,
	ShowETAcontinue1->False,
	ShowETAexpand1->False,
	ShowETAcontinue2->False,
	ShowETAexpand2->False,
	ShowMBintegrateETA->False,
	Analytical->False,
	Numerical->True
	};


OptChange[fun_,new___Rule]:=
  Module[{i},
     Map[If[Length[i=Position[{new},First[#]->_,1,1]]===0, 
            #,{new}[[Sequence @@ First[i]]]]&, Options[fun]] ]


MBnum[repr_, epsilon_, points_, powers_, loop_,options___Rule] := 
  Module[{option,AllRules, positions, cond,
          NoEtaRepr, NoEtaRules, NoEtaInts, NoEtaSer, num1list, 
          PowersEta, ReprEta, RulesEta, num2list,
	  AnalyticalResult,NumericalResult},

    
    option=OptChange[MBnum,options];

	 (*	 after= BarnesLemma[rep,1];
	  *)

 If[ShowMBrep /.option, WriteString[$Output,"repr=",InputForm[repr],
"\n","Length=",Length[repr],"\n\n"]]; 


    Rules=Map[fun[#,eps->0,{},{eps}]&, repr /.powers];

    AllRules = TimeConstrained[Rules /.fun->MBoptimizedRules,RulesConstrain /.option,Rules /.fun->MBrules];
    positions = Position[AllRules, {}, 1];


    If[ShowPositions /.option, Print["ETA's will be aplied on positions: ",Flatten[positions]]];
    If[ShowRules /.option, Print["List of all rules:\n",AllRules]];	


(*---------------------*)    
(*----- 1) NO ETA -----*)
(*---------------------*)
    If[ShowComments /.option, Print["1. Calculating 'no eta' parts..."]];	    

    NoEtaRepr = Delete[repr, positions] /.powers;
    If[Flatten[NoEtaRepr]==={},	
    	{If[ShowComments /.option, Print["   No 'no eta' parts found!!!"]];
	 NoEtaSer={};
	 num1list={}},
        {NoEtaRules= Delete[AllRules, positions];

	 If[ShowComments /.option, Print["   Running MBcontinue..."]];
    	 NoEtaInts = 
      	    MapThread[MBcontinue[#1,eps->0,#2,Verbose->False]&, {NoEtaRepr,NoEtaRules}];
         If[ShowMBcontinue /.option, Print[NoEtaInts]];


	 If[ShowComments /.option, Print["   Running MBexpand..."]];
    	 NoEtaSer = 
	    Map[MBmerge[MBexpand[{#},Exp[loop*eps*EulerGamma],{eps,0,epsilon}]]&, NoEtaInts];
	 NoEtaSer = MBmerge[NoEtaSer];
	 If[ShowMBexpand /.option, Print[NoEtaSer]];
    	
	 If[Numerical/.option,
	 If[ShowComments /.option, Print["   Running MBintegrate..."]];
         num1list = Map[MBintegrate[{#},points,Verbose->False]&, DeleteCases[NoEtaSer,{}]];
	 num1list = First[#]& /@ num1list;
	    If[ShowMBintegrate /.option, Print[num1list]]] 
	}];	
    

(*-----------------------*)    
(*----- 2) WITH ETA -----*)
(*-----------------------*)

Label["eta"];

    If[ShowComments /.option, Print["2. Calculating 'eta' parts..."]];
    If[positions==={},
	{If[ShowComments /.option, Print["   No 'eta' parts found!!!"]];
	 EtaExp={};
	 num2list={}},
	{Do[
		PowersEta = 
		   ReplacePart[powers, powers[[i]] /.Rule[n_,x_]->Rule[n,x+eta], i];
	        If[ShowComments /.option, Print["   Using eta nr: ",i,
						"\n    ",PowersEta]];

      		ReprEta = Map[Part[repr,#] /.PowersEta &, Flatten[positions]];
                RulesEta = Map[MBoptimizedRules[#,eps->0,{},{eps,eta}]&, ReprEta];

		If[ShowComments /.option, Print["   Calculating eta rules on nr: ",i]];
		If[ShowETArules /.option, Print["List of eta rules:\n",RulesEta]];


		If[ShowComments /.option, Print["   Running MBcontinue1 for eta nr: ",i]];
        	EtaInts = 
            		MapThread[MBcontinue[#1,eta->0,#2,Verbose->False]&, {ReprEta,RulesEta}];

		cond=Union[Map[MemberQ[#, MBint[___]]&,EtaInts]]//First;	

		If[cond, 
		{Break[]},
		{If[ShowComments /.option, Print["   Couldn't find rules, checking next eta..."]]}];

		If[i==Length[powers],			
			{If[ShowComments /.option, 
				Print["   Couldn't find rules for etas!!!"]];
			 Goto["end"]}],

	{i,1,Length[powers]}];


	If[ShowETAcontinue1 /.option, Print[EtaInts]];

	If[ShowComments /.option, Print["   Running MBexpand1 for eta..."]];
	EtaSer = MBexpand[#,1,{eta,0,0}]& /@ EtaInts;
	If[ShowETAexpand1 /.option, Print[EtaSer]];

	If[ShowComments /.option, Print["   Running MBcontinue2 for eta..."]];
        EtaCont = EtaSer /.MBint[integrand_,rules_]:>
                           MBcontinue[integrand,eps->0,rules,Verbose->False];
	If[ShowETAcontinue2 /.option, Print[EtaCont]];

	If[ShowComments /.option, Print["   Running MBexpand2 for eta..."]];
	EtaExp = 
      	  MBmerge[MBexpand[#,Exp[loop*eps*EulerGamma],{eps,0,epsilon}]]& /@ EtaCont;
	EtaExp = MBmerge[EtaExp];
	If[ShowETAexpand2 /.option, Print[EtaExp]];

	If[Numerical/.option,
        If[ShowComments /.option, Print["   Running MBintegrate eta..."]];
        num2list = Map[MBintegrate[{#},points,Verbose->False]&, EtaExp];
	num2list = First[#]& /@ num2list;
    	If[ShowMBintegrateETA /.option, Print[num2list]]] 
     }];

Label["end"];


AnalyticalResult = Join[NoEtaSer,EtaExp];
If[Numerical/.option,NumericalResult = Plus @@ Join[num1list,num2list]];


If[((Analytical /.option) && (Numerical /.option)),
	{result = Join[{NumericalResult},{AnalyticalResult}]},
	{If[(Analytical /.option),
		result = AnalyticalResult,
		result = NumericalResult];
	}];



Return[result]];

(* 3.06,  ShowMBrep option added
  18.06,  Numerical->False improved, adding If
*)




