(* ::Package:: *)

(*********************************************************************************************)
(*            analytic_eq.m  -  (c) 2015 Isaac Amidror   (Version 1.0,  14.6.2015)           *)
(*********************************************************************************************)
(* An interactive Mathematica(R) application accompanying the paper "Analytic derivation of  *) 
(* the (k/j)-order modulation envelopes in the sum of two mistuned (co)sine waves" by Isaac  *)
(* Amidror, 2015. This application allows one to find analytically the equation of any of    *)
(* the envelope curves that appear in the (k/j)-order modulation between two mistuned cosine *)
(* waves having frequencies f1 and f2 = (k/j)*f1 + delta. The specified envelope curve is    *)
(* plotted in red, and its analytic equation y = en(x) is displayed below the plot.          *)
(* See the detailed User Guide in the enclosed file "analytic_guide.pdf".                    *)
(*********************************************************************************************)

epsilon = 0.0000001;				(* used for comparing real numbers *)

k = 2;
j = 1;

envnumber = 1;
envmax = 2*Max[k,j];

range = 10;

f1 = 5;
delta = 0.1;
f2 = (k/j)*f1 + delta;

ratio = 2.0;
oldratio = 1000;	(* undefined *)

Off[Solve::ratnz];		(* switch off the warning message "ratnz" of the instruction Solve *)
(* Explanation: Solve, like any symbolic operation, is meant to be used with exact integer numbers
   (where roundoff errors are not an issue). Thus, whenever it gets a numeric coefficient such as
   5.25, it first converts it to (525/100), and issues the following warning message:
   "Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was 
   obtained by solving a corresponding exact system and numericizing the result."
   We prefer to switch off such annoying warnings, and not to let them pop up. *)

AbortProtect[Manipulate[
	k = 1000; j = 1000;			(* k and j are not defined *)
	If[Abs[ratio - 1.0]<epsilon, k = 1; j = 1];
	If[Abs[ratio - 2.0/1]<epsilon, k = 2; j = 1];
	If[Abs[ratio - 3.0/1]<epsilon, k = 3; j = 1];
	If[Abs[ratio - 3.0/2]<epsilon, k = 3; j = 2];
	If[Abs[ratio - 4.0/1]<epsilon, k = 4; j = 1];
	If[Abs[ratio - 4.0/3]<epsilon, k = 4; j = 3];
	If[Abs[ratio - 5.0/1]<epsilon, k = 5; j = 1];
	If[Abs[ratio - 5.0/2]<epsilon, k = 5; j = 2];
	If[Abs[ratio - 5.0/3]<epsilon, k = 5; j = 3];

	envmax = 2*Max[k,j];
	If[envnumber>envmax, envnumber=1];

	If [k >= j,
		a := Exp[2*Pi*I*delta*x];
		b := 1;
		c := -2*y;
		d := 1;
		e := Exp[-2*Pi*I*delta*x],
	(* else *)
		a := 1;
		b := Exp[2*Pi*I*delta*x];
		c := -2*y;
		d := Exp[-2*Pi*I*delta*x];
		e := 1
	];

	f2 = (k/j)*f1 + delta;

	If[ratio!=oldratio,	(* no need to solve again if ratio did not change *)
		oldratio = ratio;
		solutions := Simplify[Quiet[Solve[FullSimplify[
				Discriminant[a*t^(2*k) + b*t^(k+j) + c*t^k + d*t^(k-j) + e, t]]== 0, y]]]
	];

	penv = ListPlot[
		Flatten[
			Table[{x, y} /. {Part[solutions, envnumber]},
				{x, -range, range, ControlActive[0.5,0.1]}	(* use first the coarser step 0.5, 
																and then the finer one 0.1 *)
			], 1
		],
		Joined->True,
		PlotStyle->{Thickness[0.0035],Red},
		AxesLabel->{Style["x",Italic,12],""},
		PlotRange->{-2.1,2.1}
	];

	psum = Plot[Cos[2*Pi*f2*x] + Cos[2*Pi*f1*x],		(* plot the continuous signal s(x) *)
		{x,-range,range},
		PlotStyle->{Thickness[0.001]},
		PlotLabel->"The continuous signal \!\(\*
StyleBox[\"s\",\nFontSlant->\"Italic\"]\)(\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)) = cos(2\[Pi]\!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\)\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)) + cos(2\[Pi]\!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(2\)]\)\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)) with \!\(\*SubscriptBox[\(f\), \(1\)]\) = 5, \!\(\*SubscriptBox[\(f\), \(2\)]\) = (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))\!\(\*SubscriptBox[\(f\), \(1\)]\) + 0.1
(sum of two mistuned cosines with (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))-order modulation) and its analytically derived envelopes.
The env_n-th envelope is plotted in red; its equation \!\(\*
StyleBox[\"y\",\nFontSlant->\"Italic\"]\) = en(\!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\)) is printed below the plot.",
		AxesLabel->{Style["x",Italic,12],""},
		PlotRange->{-2.1,2.1},
		PlotPoints->751
	];

	sol1 = Quiet[Part[solutions, envnumber]];		(* pop *)

	sol = Part[sol1,1];

(*	sol = Part[Part[solutions, envnumber],1];		(* pop *)  *)

	p1 = Show[psum,penv,ImageSize->500];

	Column[{p1,Style[sol,FontFamily->"Times",10]},Center]

	,{{ratio, 2.0/1, Style["k/j",FontFamily->"Times",Italic,14]}, 
		{1.0/1->"1/1", 2.0/1->"2/1", 3.0/1->"3/1", 3.0/2->"3/2", 4.0/1->"4/1", 4.0/3->"4/3", 
		5.0/1->"5/1", 5.0/2->"5/2", 5.0/3->"5/3"}, Setter},
	{{envnumber, 1, "env_n"},Table[nn,{nn,1,2*Max[k,j],1}],Setter},
	{{j,1,8},ControlType->None},	(* dummy, but necessary! *)
	TrackedSymbols->All,
	SaveDefinitions->True,		 (* required for generating a CDF file *)
	ContinuousAction->False		(* to avoid the 5-second limit after which "$Aborted" is displayed *)
]]

