;

Print[ Text[Style["\n The Mathematica code can be obtained from the author upon request.  ",Darker[ Blue],Italic,28]]];

(* This code is provided by the author, I. Green, to complement the subject paper. You may use the code freely. *)
(* Exercise any or all of the cases listed in Table 1 in the said paper. This code will reproduce all of the results in there. *)
(* Missing from this code is the implementation of the spectral moments as found in Ref.[10] because they produce results that are grossly in error. *)
(* The errors are caused by the approximated power spectrum, sometimes dubbed as "the continuous power spectrum density," as derive in Ref. [11]. 
That spectrum had been used in Ref. [10] and, since then, in many other papers. The errors sometimes approach 100%!  
These errors put the "approximated/continuous power spectrum," and the so-called "power law," in question.  *)

(* NOTE: dmp and dmp24 are proprietary print statements. Use Mathematica's "Print[]" statement instead. *)

dothis=1; (* switch - do analytical integrations of the spatial signal using their derivatives *)
DD=1.5;g=1.5;G=1.0; (* fractal parameters;  *)
dmpc24[{DD,g,G},Darker[Red]];

LL=100; (* The spectral moments should be independent of LL; Indeed, they don't much when LL is large, and the differentiation method is used; It doesn't at all in EXACT soltion! *)
L0=LL*0/5; (* arbitrary offset *)

(* change switch q and see how units affect the results *)
q=1*1; (* q-units switch: If =1 then [Hz]; If =0 then [rad/s] *)
dmpc24[q,Darker[Red]];

n1=0;
R[x_]:=(G^(DD-1)Sum[Cos[(2*Pi)^q*g^i*(x)]/g^(i*(2-DD)),{i,n1,n2}] // Expand)
m[k_] :=((G^(2*DD-2)/2)   *Sum[(2*Pi)^(q*k)/g^(i*(4-2*DD-k)),{i,n1,n2}] // Expand)     (* These are the exact spectral moments! switched by q  *)

color={Blue,Red,Magenta,Brown,Orange,Pink,Purple,Black};
N2={17,26,34}/ 1 // Floor; (* cases of n2 you want to run; These values are taken from Table 1 *)
noplots=Length[N2]; (* number of cases to run; below see: j=1 to noplots *)
Do[Print["\n"];
n2=N2[[j]]; (* this parameter value DOES change spectral moments m2 and m4, but not m0 !! *)
dmpc24[n2,color[[j]]];
chopplotfactor=10;  (* make chopplotfactor=1 to see the whole signal range; if greater than zero, the range is chopped in the plot only *)
gr[j]=Plot[R[x+L0],{x,0,LL/chopplotfactor}, GridLines->Automatic, PlotStyle->{color[[j]]}];Print[gr[j]];

Print[ Text[Style["\n The EXACT Spectral Moments according to the IG paper (it is implied that LL->Infinity); \n These are entirely analytical, so there is NO CPU time consumed, at all! ",ColorData["Atoms"]["Ba"],Italic,28]]];
dmpc24[sm0=m[0],ColorData["Atoms"]["Ba"]];
dmpc24[sm2=m[2],ColorData["Atoms"]["Ba"]];
dmpc24[sm4=m[4], ColorData["Atoms"]["Ba"]];

If[dothis==1,
tm=DateString[];t0=AbsoluteTime[tm];
avg=(1/LL)*Integrate[R[x+L0] // Expand,{x,0,LL}]; dmp[avg];
m0=(1/LL)*Integrate[(R[x+L0]-avg)^2 // Expand,{x,0,LL}];
m2=(1/LL)*Integrate[D[R[x+L0],x]^2   // Expand,{x,0,LL}];
m4=(1/LL)*Integrate[D[R[x+L0],{x,2}]^2  // Expand  ,{x,0,LL}];
Print[ Text[Style["\n Spectral Moments by Differentiation (the results are exact for this signal length, LL);\n This method does consume CPU time! The larger n2, the larger the CPU time!", Blue,Italic,28]]];
dmpc24[{m0,m2,m4}, Darker[Blue]];
tm=DateString[];t1=AbsoluteTime[tm];dmp[cpu[seconds]=t1-t0];

Print["\n"];
Print[ Text[Style["\n Relative difference between differentiation method and exact solution \n(the larger LL, the smaller the difference):  ",ColorData["Atoms"]["Ba"],Italic,28]]];
prt0["reldiff0=Abs[sm0-m0]/sm0"];
prt0["reldiff2=Abs[sm2-m2]/sm2"];
prt0["reldiff4=Abs[sm4-m4]/sm4"]; 
Print["\n "]];

,{j,1,noplots}]

Print[ Text[Style["\n Fractal plots of all cases above ", Blue,Italic,28]]];

lst=Table[gr[j],{j,noplots,1,-1}];  (* after the do-loop *)
Show[lst]