(* ::Package:: *)

(*********************************************************************************************)
(*           analytic_env.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 dynamically observe beating effects that    *)
(* occur between two mistuned cosine waves (or between two mistuned sine waves) having       *)
(* frequencies f1 and f2 = (k/j)*f1 + delta. The button "env" allows to select the desired   *)
(* envelope type: none, cosinusoidal approximations (braid curves), or analytically derived  *)
(* envelope curves, as explained in the paper.                                               *)
(* The value of the frequency f1 can be freely swept between 0.5 and 5, in order to show how *)
(* f1 influences the signal's beats and their envelopes.                                     *)
(* The relative phase phi can also be varied between 0 - 1, giving various possible effects  *)
(* as determined by the button "phase" (see the explanations provided for mistuned.m).       *)
(* See the detailed User Guide in the enclosed file "analytic_guide.pdf".                    *)
(*********************************************************************************************)

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

f1 = 1.0;			(* frequency of the first cosine *)
T1 = 1.0/f1;		(* period of the first cosine *)

range = 16;			(* signal range: -range...range *)

envelopecol1 = {Red,Green,Blue,Cyan,Magenta,Orange,Black,Yellow};		(* color sequence for the upper envelopes *)
envelopecol2 = {Cyan,Magenta,Orange,Black,Yellow,Red,Green,Blue};		(* color sequence for the lower envelopes *)

ratio = 2.0/1;		(* the integer ration k/j *)
del = 1.0/8;		(* the distance delta between f2 and (k/j)*f1 *)
f2 = ratio*f1+del;	(* frequency of the second cosine: can vary between 0...2.5f1 *)

env = 0;		(* 1 = plot the envelope approximations; 2 = plot the analytical envelopes; 0 = don't plot the envelopes *)
highqualityplot = True;	(* False = plot curves using lower PlotPoints, which is quicker but lossy *)

oldf1 = f1;
oldphi = 0;

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[
	If[highqualityplot, pointdistance = 0.25, pointdistance = 1];

	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 - 4.0/1]<epsilon, k = 4; j = 1];
	If[Abs[ratio - 3.0/2]<epsilon, k = 3; j = 2];

	f2 = ratio*f1 + del;
	T2 = If[f2==0, 1000000, 1/f2];		(* period of the second cosine *)

	fenvbc = If[k>j, del*j/k, If[k==j, del/2, del]];		
													(* frequency of the (k/j)-bc-envelope *)
	Tenvbc = If[fenvbc==0, 1000000, 1/fenvbc];	(* period of the (k/j)-bc-envelope *)
	fenv = del*j/(k+j);					(* frequency of the (k/j)-envelope; see mistuned.m *)
	Tenv = If[fenv==0, 1000000, 1/fenv];	(* period of the (k/j)-envelope; see mistuned.m *)

	Switch[phasetype,
		0,	(* rigid shift: *)
			  phi1 = f1*phi*Tenvbc;		(* multiply by Tenvbc to increase the rigid shift *)
			  phi2 = f2*phi*Tenvbc;		(* multiply by Tenvbc to increase the rigid shift *)
			  phienv = phi;				(* phienv - phase of the (k/j)-envelope; see mistuned.m *)
			  phienvbc = phi,			  (* phienvbc - phase of the (k/j)-bc-envelope *)
		1,	(* only envelopes move: *)
			  phi1 = -phi;
			  phi2 = phi;
			  phienv = phi;
			  phienvbc = phi*fenvbc/fenv,
		2,	(* only carriers move: *)
			  phi1 = j*phi;
			  phi2 = k*phi;
			  phienv = 0;
			  phienvbc = 0,
		3,	(* both carriers and envelopes move, in the same direction and speed: *)
			  phi1 = j*phi*Tenvbc - phi;	(* multiply by Tenvbc to increase the carriers' shift *)
			  phi2 = k*phi*Tenvbc + phi;	(* multiply by Tenvbc to increase the carriers' shift *)
			  phienv = 0 + phi;
			  phienvbc = 0 + phi*fenvbc/fenv,
		4,	(* both carriers and envelopes move in the same speed, but in opposite directions: *)
			  phi1 = -j*phi*Tenvbc - phi;	(* multiply by Tenvbc to increase the carriers' shift *)
			  phi2 = -k*phi*Tenvbc + phi;	(* multiply by Tenvbc to increase the carriers' shift *)
			  phienv = 0 + phi;
			  phienvbc = 0 + phi*fenvbc/fenv
		(* Note that in cases 3 and 4 the first term of phi1 and of phi2 corresponds to
			the phase of the carrier (see case 2), while the second term corresponds to
			the phase of the envelopes (see case 1). *)
	];

	(* Generate the top figure: *)

	p1a = Plot[				(* plot the continuous signal s(x) *)
			Cos[2*Pi*f2*x + 2*Pi*phi2] + Cos[-2*Pi*f1*x - 2*Pi*phi1],
		{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[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\), \!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(2\)]\) = (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))\!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\) + \[Delta]
(sum of two mistuned cosines with (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))-order modulation)",
		AxesLabel->{Style["x",Italic,12],""},
		PlotRange->{-4.01,4.01},
		PlotPoints->751			(* reducing this value will speed up the application, 
									at the expense of possibly missing some of the rapid
									oscillations of the signal s(x) being plotted here *)
	];

	n = Max[k,j];

	(* if env \[Equal] 1 plot the cosinusoidal envelope approximations (braid curves): *)

	If[env==1 && k+j<500,			(* plot the (k/j)-order envelopes only if k,j are defined *)
		penv = {};
		For[m=0, m<n, m++,		(* there are 2n interlaced envelopes (braid curves) *)
			pe1 = Plot[			(* plot the n envelopes above the x axis *)
				If[k==1 && j==1, 2*Cos[2*Pi*fenvbc*x + 2*Pi*phienvbc],
						Cos[2*Pi*fenvbc*x + 2*Pi*(phienvbc+(m*j/k))] + 1		(* else: k>j *)
				],
				{x,-range,range},
				PlotStyle->{Thickness[0.0045],envelopecol1[[m+1]]},
				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[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\), \!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(2\)]\) = (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))\!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\) + \[Delta]
(sum of two mistuned cosines with (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))-order modulation)",
				AxesLabel->{Style["x",Italic,12],""},
				PlotRange->{-4.01,4.01},
				PlotPoints->351
			];
			AppendTo[penv,pe1];
			pe2 = Plot[			(* plot the n envelopes below the x axis *)
				If[k==1 && j==1, -2*Cos[2*Pi*fenvbc*x + 2*Pi*phienvbc],
						Cos[2*Pi*fenvbc*x + 2*Pi*(phienvbc+((m-0.5)*j/k))] - 1	(* else: k>j *)
				],
				{x,-range,range},
				PlotStyle->{Thickness[0.0045], If[k==1 && j==1, envelopecol1[[m+2]], envelopecol2[[m+1]]]},
				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[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\), \!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(2\)]\) = (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))\!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\) + \[Delta]
(sum of two mistuned cosines with (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))-order modulation)",
				AxesLabel->{Style["x",Italic,12],""},
				PlotRange->{-4.01,4.01},
				PlotPoints->351
			];
			AppendTo[penv,pe2];
		];					(* end of For *)
		p1 = Show[p1a,penv],	(* penv contains the 2n envelope approximations (braid curves) *)

		(* else if env\[Equal]2 plot the analytically derived envelopes: *)

		If[env==2 && k+j<500,		(* plot the (k/j)-order envelopes only if k,j are defined *)
			If[oldf1 == f1,			(* if the slider f1 was moved use old penv, don't recalculate *)
				penv = {};
				solutions := Quiet[Solve[
					Discriminant[Exp[2*Pi*I*del*(x + phienvbc*Tenvbc)]*t^(2*k) +
						1*t^(k+j) + (-2*y)*t^k + 1*t^(k-j) + 
						Exp[-2*Pi*I*del*(x + phienvbc*Tenvbc)], t] == 0, y]];
				For[m=0, m<2*n, m++,		(* there are 2n envelopes; see Theorem *)
					pe = ListPlot[
						Flatten[
							Table[{x, y} /. {Part[solutions, m+1]},
								{x, -range, range, pointdistance}
							], 1
						],
						Joined->True,
						PlotStyle->{Thickness[0.0045],envelopecol1[[2*n - m]]},
						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[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\), \!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(2\)]\) = (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))\!\(\*SubscriptBox[
StyleBox[\"f\",\nFontSlant->\"Italic\"], \(1\)]\) + \[Delta]
(sum of two mistuned cosines with (\!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)/\!\(\*
StyleBox[\"j\",\nFontSlant->\"Italic\"]\))-order modulation)",
						AxesLabel->{Style["x",Italic,12],""},
						PlotRange->{-4.01,4.01}
					];
					AppendTo[penv,pe];
				]			(* end of For *)
			];
			p1 = Show[p1a,penv],	(* penv contains the 2n analytically derived envelopes *)
			p1 = Show[p1a]		(* else (i.e. if env\[Equal]0 or k,j undefined) do not plot envelopes *)
		]
	];
	oldf1 = f1;

	(* Generate the bottom figure: *)

	p2 = Plot[				(* plot the two original mistuned cosines *)
		{Cos[2*Pi*f1*x + 2*Pi*phi1]*If[x>12.5,0,1],
		 Cos[2*Pi*f2*x + 2*Pi*phi2]*If[x<-12.3,0,1]},
		{x,-range,range},
		PlotStyle->{{Thickness[0.001],Cyan},{Thickness[0.001],Magenta}},
		PlotLabel->"The two given mistuned cosines (Cyan and Magenta)",
		AxesLabel->{Style["x",Italic,12],""},
		PlotRange->{-4.01,4.01},
		PlotPoints->751			(* reducing this value will speed up the application, 
									at the expense of possibly missing some of the rapid
									oscillations of the signal s(x) being plotted here *)
	];

	pp = GraphicsGrid[{{p1},{p2}},ImageSize->500]

	,Row[{Item["   "], 
		Control[{{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", 
		4.0/1->"4/1", 3.0/2->"3/2"}, Setter}],
		Item["                                    "],
		Control[{{highqualityplot, True, "plot density"}, {True->"high", False->"low"}}],
	}],
	{{del, 1.0/8, Style["d",FontFamily->"Symbol",Italic,14]}, {-1.0/8->"-1/8", -1.0/16->"-1/16",
		0->"0", 1.0/16->"1/16", 1.0/8->"1/8"}, Setter},
	Row[{Item["   "], 
		Control[{{env, 0, "env"}, {0->"none", 1->"approx.", 2->"analytic"}}],
		Item["     (which envelope type is to be plotted)"]
	}],
	{{f1, 1, Subscript[Style["f",FontFamily->"Times",Italic,14],Style["1",FontFamily->"Times",8]]}, 0.5, 5, Grid[{{
				Slider[##, ImageSize->400],InputField[#, FieldSize->Tiny]
	}}, Alignment->{Center,Center}] &},
	{{phasetype, 0, "phase"}, {0->"rigid", 1->"envelope", 2->"carrier", 3->"both", 4->"divergent"}},
	{{phi, 0, Style["f",FontFamily->"Symbol",Italic,14]}, 0, 1, Grid[{{
				Slider[##, ImageSize->400],InputField[#, FieldSize->Tiny]
	}}, Alignment->{Center,Center}] &},
	TrackedSymbols:>{ratio,f1,del,phasetype,phi,env,highqualityplot},
	SaveDefinitions->True				(* required for generating a CDF file *)
]]

