User:Deciwill

// Optomechanical Theory // V 1.1 D Brooks //	optomechanical theory rewritten to work with global variables, fit-compatible functions, and a front panel. //	Since math and pi cannot be entered into panel, we will work in cyclic units. //	0 phase shift is antinode. node is at pi*0.5. // v 1.3	N Brahms //	Modified for 9-well averaging // v5		N Brahms //	Integrated with version 5 code, arbitrary well averaging, some cleanup.
 * 1) pragma rtGlobals=1		// Use modern global access method.

function omFunctionList // PANELS and INITIALIZATION //	makeOmPanel // 	omInit // CALCULATE OM PARAMETERS //	omRecalc				//	calcAParameters //	myprotofunc				//	avgWell //	getContrast				//	getDN //	getKOpt					//	BroadenSlope //	getStaticK				//	getStaticOsc //	getOscFreq				//	getGQed //	getZHo					//	getLambDicke //	getGranularity			//	getHeating //	gOpt // FITTING FUNCTIONS //	FreqVsDelta_fit			//	FreqVsPhi_fit //	fano						//	logFano //	BreitWigner				//	BreitWigner2 //	omFitMaster				//	SSresponse end

Menu "Macros" "Open Optomechanics Panel", makeOmPanel end

Window makeOmPanel: Panel DoWindow /k omPanel NewPanel /w=(370,320,370+520,320+230) /N=omPanel /k=1 as "Optomechanics Panel" //Optomechanic variables variable/g OptoMechInit		//Initializes OptoMechanic variables only first time panel is called if(!OptoMechInit) omInit("") endif //Other parameters //variable/g sLa0, sLa1p, sLa1n, sLa2p,sLa2n,sLn 		// linear sensitivity for five well weightings make/o/t/n=(numpnts(omA)) omAList omAList = "A["+num2str(p)+"] = "+num2str(omA[p]) // variable i	// for (i=0; i<numpnts(omA); i+=1) //	omAList[p] = num2str(omA[p]); // endfor //column 1 variable c=10;	variable r=5; button resetparams,	pos={c,r},	size={150,30},	proc=omInit,			title="Click to Reset\r Parameters to Default";	r+=40 setvariable omG0,	pos={c,r},	size={150,15},	proc=omRecalc,		title="g_0 [Hz] ",			value=omG0,			bodywidth=70,	format="%4.3g",	disable=2;	r+=20 setvariable omK,		pos={c,r},	size={150,15},	proc=omRecalc,		title="kappa [Hz] ",		value=omK,			bodywidth=70,	format="%4.3g",	disable=2;	r+=30 setvariable omDCa,	pos={c,r},	size={150,15},	proc=omRecalc,		title="Delta_ca [Hz] ",		value=omDCa,		bodywidth=70,	format="%4.3g";	r+=20 setvariable omDelta,	pos={c,r},	size={150,15},	proc=omRecalc,		title="Delta [Hz] ",		value=omDelta,		bodywidth=70,	format="%4.3g";	r+=20 setvariable omNBar,	pos={c,r},	size={150,15},	proc=omRecalc,		title="Num photons ",		value=omNBar,		bodywidth=70;	r+=20 setvariable omNPhonon,	pos={c,r},	size={150,15},	proc=omRecalc,	title="Num phonons ",		value=omNPhonon,	bodywidth=70;	r+=20 setvariable omNA,	pos={c,r},	size={150,15},	proc=omRecalc,		title="Num atoms ",		value=omNA,			bodywidth=70;	r+=20 setvariable omf0,pos={c,r},	size={150,15},	proc=omRecalc,		title="f 0 [Hz] ",			value=omf0,	bodywidth=70,	format="%4.3g";	r+=20

//column 2 c+=170; r=5; setvariable omPhi,	pos={c,r},	size={150,15},	proc=omRecalc,		title="Phi [pi] ",			value=omPhi,		bodywidth=70;	r+=20 setvariable omNoise,	pos={c,r},	size={150,15},						title="Cavity noise [Hz] ", 	value=omDeltaNoise,	bodywidth=70,	format="%4.3g";	r+=20; setvariable omASize,	pos={c,r},	size={150,15},	proc=omRecalc,	title="Dist. width [a]",		value=omASize,		bodywidth=70;	r+=20 setvariable omNumA,	pos={c,r},	size={150,15},	proc=omRecalc,		title="# avg. coeffs. ",		value=omNumA,		bodywidth=70;	r+=30 valdisplay omContrast,pos={c,r},	size={140,15},						title="Contrast: ",			value=omContrast,	bodywidth=50,	format="%4.3f";	r+=20 valdisplay omDN,		pos={c,r},	size={140,15},						title="Delta_N [MHz]: ",	value=omDN,			bodywidth=50,	format="%4.1f";	r+=25 listbox omAList,		pos={c+20,r},	size={130,72},					listWave = omAList; // valdisplay omA0,pos={c+5,r},size={85,15},title="a0 ",value=omA0;r+=20 // valdisplay omA1,pos={c+5,r},size={85,15},title="a1 ",value=omA1;r+=20 // valdisplay omA2,pos={c+5,r},size={85,15},title="a2 ",value=omA2;r+=20 //column 3 c+=180; r=5 titlebox omFreqTitle,	pos={c,r},	size={150,15},	frame=0,		title="Oscillation frequencies (Hz):";									r+=20 valdisplay omOmegaOsc,pos={c,r}, size={150,15},				title="Total ", 				value=omOmegaOsc,	bodywidth=60,	format="%4.3g";	r+=20 valdisplay omOmegaStatic,pos={c,r}, size={150,15},			title="Static ", 				value=omOmegaStatic,bodywidth=60,	format="%4.3g";r+=20 valdisplay omOmegaShift,pos={c,r}, size={150,15},			title="Shift ", 					value=omOmegaShift,	bodywidth=60,	format="%4.3g";	r+=30 valdisplay omGQed,	pos={c,r},	size={150,15},				title="cQED coupling (Hz) ",	value=omGQED,		bodywidth=60,	format="%4.3g";	r+=20 valdisplay omGOm,	pos={c,r},	size={150,15},				title="OM coupling (Hz) ",		value=omGOm,		bodywidth=60,	format="%4.3g";	r+=20 valdisplay omZSize,	pos={c,r},	size={150,15},				title="Z size (nm) ",			value=omZSize,		bodywidth=60,	format="%4.1f";	r+=20 valdisplay omEpsilon,	pos={c,r},	size={150,15},			title="Granularity param ",		value=omEpsilon,		bodywidth=60,	format="%4.3f";	r+=20 valdisplay omHeating,	pos={c,r},	size={150,15},			title="Heating (quanta/s) ",		value=omHeating,		bodywidth=60,	format="%4.3g";	r+=20

omRecalc("",0,"","") end

//return to omFunctionList //---

//--- function omInit(ctrlname) : ButtonControl String ctrlName variable/g omG0 = 13.1e6				// per-atom on resonance coupling strength variable/g omK = 1.82e6				// cavity linewidth via ringdown measurement variable/g omPhi	= 1/4				// probe-trap phase shift at atom well, 0 is antinode variable/g omNBar = 3.5				// mean photon number, averaged over the sidelock bandwidth variable/g omNA	= 4000				// number of atoms variable/g omNPhonon = 0.15			// number of thermal COM phonons variable/g omDCa = -40e9				// detuning of bare cavity resonance from atomic resonance (red is negative) variable/g omMass = 1.44466881e-25	// per-atom mass of Rb, m=87 amu variable/g omKP	= 2*pi/(780e-9)		// probe wavenumber variable/g omKT	= 2*pi/(850e-9)		// trap wavenumber variable/g omf0	 = 57.7e3		// trap frequency variable/g omfs	= 140e3				// shifted frequency variable/g omDelta = -2.05e6			// detuning of probe from averaged with-atom cavity resonance (red is negative) variable/g omASize = .88				// Size of cloud width in units of trap wells (425 nm) variable/g omNumA = 3				// Number of wells to include in well-weighting (must be integer > 1) make/o/n=(omNumA) omA				// Site averaging coefficients make/o/n=(2*omNumA-1) omL			// Linear sensitivity coefficients //variable/g omA0,omA1,omA2,omA3,omA4 variable/g omDPhi = 0.27	/pi			// phi spacing for multi-well averaging; There are 11.64 wells per phase of pi variable/g omDeltaNoise = .6e6			// Standard Deviation of the gaussian technical noise that broadens cavity linewidths and slopes. variable/g omContrast variable/g omOmegaOsc variable/g omOmegaStatic variable/g omOmegaShift variable/g omGQed					// Effective cavity QED coupling parameter, = sqrt(N)*g0^2/dca variable/g omGOm					// Effective optomechanical coupling parameter, = sqrt(N)*g0^2/dca*kp*zho variable/g omZSize					// Size within well (nm) variable/g omEpsilon					// Granularity parameter variable/g omHeating					// Per-atom heating rate (quanta/s) variable/g omDN						// Delta N (MHz) omRecalc("",0,"","") nvar OptoMechInit OptoMechInit=1 end

//return to omFunctionList // //	function omRecalc(ctrlName, varNum,varStr, varName) //		recalculates omContrast and linear sensitivity parameters when a variable is changed on the front panel //		can also be called from functions, useful during fit routines // function omRecalc(ctrlName, varNum,varStr, varName) :SetVariableControl String ctrlName Variable varNum	// value of variable as number String varStr		// value of variable as string String varName	// name of variable calcAParameters nvar omContrast, omASize, omPhi omContrast=getContrast nvar omOmegaOsc, omOmegaStatic, omOmegaShift, omDelta,omNA,omf0 nvar omGQed, omGOm, omZSize, omEpsilon, omHeating, omKP, omNPhonon, omDN omOmegaOsc = getOscfreq(omDelta,omNA,omf0,omASize,omPhi,initialized=1) omOmegaStatic = getStaticOsc(omPhi) omOmegaShift=Sign(omOmegaOsc^2-omOmegaStatic^2)*(Abs(omOmegaOsc^2-omOmegaStatic^2))^.5 omGQed = getGQed; omGOm = omKP*getZHO*omGQed*sqrt(omNA); //	omLD = getLambDicke; omZSize = (2*omNPhonon+1)*getZHO*1e9; omDN = avgWell(getDN)/1e6; omEpsilon = getGranularity; omHeating = avgWell(getHeating);

DoUpdate /W=omPanel end

//return to omFunctionList // // function calcAParameters // function calcAParameters nvar omASize,omNumA,omPhi,omDPhi // Need to average over at least one well if (mod(omNumA,1)!=0 || omNumA < 1) DoAlert 0, "Number of wells must be an integer >= 1" wave omA omNumA = numpnts(omA) return -1 endif // A size must be positive if (omASize < 0) DoAlert 0, "A size must be positive" wave omA;							omA = NaN; make/o/t/n=(numpnts(omA)) omAList;	omAList = "A["+num2str(p)+"] = NaN"; return -2 endif make/o/n=(omNumA) omA // Calculate the well-weighting parameters if (omASize==0) omA = 0;	omA[0] = 1; elseif (omNumA>1) variable iA, aTotal = 0 for(iA=0; iA<omNumA-1; iA+=1) omA[iA] = .5*Erf((iA+.5)/omASize/sqrt(2)) - aTotal; aTotal += omA[iA]; endfor omA[iA] = .5*Erfc((iA-.5)/omASize/sqrt(2));	// Last wells omA[0] *= 2;								// Correct for only having one center well else omA[0] = 1; endif // Make labels for the list box make/o/t/n=(numpnts(omA)) omAList omAList = "A["+num2str(p)+"] = "+num2str(omA[p]) // Calculate linear sensitivity make/o/n=(2*omNumA-1) omL omL[omNumA-1] = omA[0]*sin(2*pi*omPhi)		// Handle center well for(iA=1; iA<numpnts(omA); iA+=1) omL[omNumA-iA-1] = omA[iA]*sin(2*pi*(omPhi-iA*omDPhi)) omL[omNumA+iA-1] = omA[iA]*sin(2*pi*(omPhi+iA*omDPhi)) endfor end

//return to omFunctionList // // function myprotofunc(localPhi) //	serves as a template for functions which are to be passed to another function using FUNCREF //	takes a phase 0 to 0.5 as an input //	returns some calculation based on that omPhi // function myprotofunc(localPhi) variable localPhi nvar omDPhi return localPhi+omDPhi end

// // function avgWell(fin,[localPhi]) // 	Averages input function, fin, over omNumA wells weighted by atom density in those wells //	if a omPhi isn't specified, uses the global omPhi value // function avgWell(fin,[phi]) FUNCREF myprotofunc fin variable phi nvar omDPhi,omPhi wave omA

//Default Parameters if(ParamisDefault(phi)) phi=omPhi endif variable avg = omA[0]*fin(phi)			// Center well variable iA	for (iA=1; iA<numpnts(omA); iA+=1)	// All other wells avg += omA[iA]*(fin(phi+iA*omDPhi)+fin(phi-iA*omDPhi)) endfor return avg end //return to omFunctionList

// // getContrast returns the 5 well averaged omContrast // function getContrast variable DNmax, DNmin DNmax = avgWell(getDN,phi=0) DNmin = avgWell(getDN,phi=0.5)

return (DNmax-DNmin)/(DNmax+DNmin) end // // getDN calculates the atom-cavity shifted resonance for the atoms at phase localphi // function getDN(phi) variable phi variable DNphipos nvar omNA,omG0,omDCa,hbar,omKP,omMass,omf0,omNPhonon DNphipos=(omNA*omG0^2/omDCa)*.5*( 1 + exp(-getLambDicke^2) * cos(2*pi*phi) ) return DNphipos end //return to omFunctionList

// // getKOpt calculates the oscillation frequency of the atoms located at local omPhi // function getKOpt(localPhi) variable localPhi variable statick, dynamickprefactor, dynamick nvar 	omDelta,omNA,omf0,omASize,omPhi nvar hbar, omG0, omKP, omDCa, omNBar, omDPhi,omK,omMass,omDeltaNoise

statick = -2 * hbar * (2*pi*omG0)^2 * omKP^2 * omNBar / (2*pi*omDCa)*cos(2*pi*LocalPhi) dynamickprefactor = 2 * hbar * omKP^2 * omNA * (2*pi*getGQed)^2 * omNBar* sin(2*pi*LocalPhi)^2 dynamick=dynamickprefactor*Integrate1d(BroadenSlope,(-5*omDeltaNoise),(5*omDeltaNoise)) return (statick+dynamick) end

// //BroadenSlope is used to calculate noise-broaden slope to be used for dynamic shifts // function BroadenSlope(inx) variable inx nvar omDelta, omK, omDeltaNoise variable slopefunc=((2*pi*(omDelta+inx)) / ((2*pi*omK)^2 + (2*pi*(omDelta+inx))^2)	return 1/sqrt(2*pi)/(omDeltaNoise)*Exp(-(inx^2)/(2*omDeltaNoise^2))*slopefunc end //return to omFunctionList

// // getStaticK calculates the static shifts only // function getStaticK(localPhi) variable localPhi variable statick, dynamick nvar 	omDelta,omNA,omf0,omASize,omPhi nvar hbar, omG0, omKP, omDCa, omNBar, omDPhi,omK,omMass statick = -2 * hbar * (2*pi*omG0)^2 * omKP^2 * omNBar / (2*pi*omDCa)*cos(2*pi*LocalPhi) return statick end

// // getStaticOsc calculates the oscillation frequency from static shifts only // function getStaticOsc(localPhi) variable localPhi nvar omf0,omMass variable f=sqrt((2*pi*omf0)^2+(avgWell(getStaticK,phi=localPhi))/omMass)/2/pi return f end //return to omFunctionList

//- //function getOscFreq(LocDelta,LocNa,LocWt,LocAsize,LocPhi,[initialized]) //	Returns the cyclic frequency of the oscillator, //	as a function of omDelta, # of atoms, trap frequency, distribution of atoms, and omPhi. // 	Returns 5-well sensitivity averaged results. // //	To be used by fit function, it needs to take local variables, which are then written over the global variables //	initialized = 0 by default. use initialized =1 when called from omRecalc //- function getOscFreq(LocDelta,LocNa,LocWt,LocAsize,LocPhi,[initialized]) variable LocDelta,LocNa,LocWt,LocAsize,Locphi variable initialized //Default Parameters if(ParamisDefault(initialized)) initialized=0 endif nvar 	omDelta,omNA,omf0,omASize,omPhi omDelta=LocDelta; omNA=LocNa; omf0=LocWt; omASize=LocAsize; omPhi=LocPhi; if(! initialized) omRecalc("",0,"","") endif nvar omMass,omf0

variable f=sqrt((2*pi*omf0)^2+(avgWell(getKOpt))/omMass)/2/pi return f end //return to omFunctionList

// // function getGQed //	Calculates the single-atom cavity QED coupling parameter in cyclic frequency //	units. // //	gQED = g0^2 / Dca //

function getGQed nvar omNA, omG0, omDca return omG0^2/omDca end

// // function getZHo //	Calculates the (single-atom) harmonic oscillator length // //	zHO = sqrt(h / 2 m f0) //

function getZHo nvar hbar, omMass, omf0 return sqrt(hbar/omMass/4/pi/omf0) end //return to omFunctionList

// // function getLambDicke //	Calculates the Lamb-Dicke parameter // //	LD = sqrt(2*nth+1)*zHO*kt //

function getLambDicke nvar omNPhonon, omKT return sqrt(2*omNPhonon+1)*getZHo*omKT end

// // function getGranularity //	Calculates the (many-atom) granularity parameter using global //	values of optomechanics parameters. // //	Uses //		epsilon = sqrt(N) F_dip z_HO / hbar kappa //	where //		F_dip = hbar kp g0^2 / dca // // c.f. Murch thesis //-

function getGranularity nvar omKP, omK, omNA return omKP*getZHO*sqrt(omNA)*getGQed / omK end //return to omFunctionList

//- // function getHeating //	Returns the per-atom heating rate, in quanta / s // //	Assumes heating only from Snn^(-) // //	Heating is therefore 2*pi*nbar*kappa*epsilon^2/N //--

function getHeating(phi) variable phi nvar omK, omNBar, omNA variable ld 	ld = sqrt( sin(2*pi*phi)^2 + getLambDicke^2*cos(2*pi*phi)^2 ) return 4*pi^2*omNBar*omK*getGranularity^2/omNA*ld	// omK is in cyclic units end //--- //function gOpt calculates OM damping rate function gOpt(freq) variable freq nvar omK, omf0 return 2 * omK * (omf0^2 - freq^2) / W1(freq) end //function W0 = kappa^2 + delta^2 function W0 nvar omK, omDelta end

// function W1 = (kappa^2 + delta^2 - freq^2 ) function W1(freq) variable freq nvar omK, omDelta return omK^2 + omDelta^2 - freq^2 end // function ksf not finished function ksf(freq) variable freq nvar omfs return (omfs^2 - freq^2)* W0/W1(freq) end //--- //function FreqVsDelta_fit(w,t) : fitfunc //	Fitting Function that calls GetOscFreq //	The x wave should be a list of deltas. //	Uses cyclic frequencies for parameters //-

function FreqVsDelta_fit(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 6 //CurveFitDialog/ w[0] = omDelta_scaling //CurveFitDialog/ w[1] = omDelta_offset //CurveFitDialog/ w[2] = AtomNum //CurveFitDialog/ w[3] = Unshifted_TrapFreq //CurveFitDialog/ w[4] = AtomSize //CurveFitDialog/ w[5] = omPhi return GetOscFreq(w[0]*t+w[1],w[2],w[3],w[4],w[5])

end //return to omFunctionList

//--- //function FreqVsPhi_fit(w,t) : fitfunc //	Fitting Function that calls GetOscFreq //	The x wave should be a list of omPhi's //	Uses cyclic frequencies for parameters //-

function FreqVsPhi_fit(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = omDelta //CurveFitDialog/ w[1] = AtomNum //CurveFitDialog/ w[2] = Unshifted_TrapFreq //CurveFitDialog/ w[3] = AtomSize //CurveFitDialog/ w[4] = omPhi_Offset return GetOscFreq(w[0],w[1],w[2],w[3],t+w[4])

end

//--- // FITTING FUNCTIONS //---

//--- //function omFitMaster(c,f,mode,[doLog]) // //	Master function for all OM response fits. Has variable //	functionality depending on choice of "mode" flag (see function //	comments). For modes 0 and 1, returns a power spectral density //	as a function of frequency "f". For mode 2 and 3, returns the real //	and imaginary parts, respectively, of the response to transmitted //	AM drive. // //	Note: This function is NOT a fitting function! //	Use a derived function. //--- function omFitMaster(c,f,[mode,doLog]) wave c			// c is the coefficients wave. Note that all frequencies must have identical units //	(i.e. rad/s, kHz, etc.) // WARNING: In order to preserve the functionality of the various functions that //	use omFitMaster, the order of these coefficients should NOT be changed //	without changing all dependent functions. Consider adding new coefficients //	at the END of the existing list. // c[0] -- Natural frequency // c[1] -- Shifted frequency // c[2] -- Natural damping rate (frequency units) // c[3] -- Cavity half-linewidth (frequency units) // c[4] -- Cavity detuning (frequency units) // c[5] -- Shot noise level (used only for modes 0 & 1) // c[6] -- Maximum technical noise level //			For mode 0 -- Units of noise level //			For modes 2-3 -- AM drive level // c[7] -- Minimum technical noise level, (noise level, used only for mode 0) // c[8] -- Noise peak quadrature angle (radians, used only for mode 0) // c[9] -- Detection quadrature angle (used only for modes 0 & 1) // c[10] -- Detection efficiency (used only for modes 0 & 1) variable f			// f is the response frequency variable mode	// Controls operation of omFitMaster // Possible values of mode (default is 0): // 0 - Normal operation, PSD driven by shot noise and constant technical noise // 1 - PSD driven by shot noise and technical noise in waves "techAm" and "techFm" //		Due to subtleties of technical noise, this mode can not be used for quadratures other than 0 or pi/2. // 2 - Network analysis, real response, driven by AM noise before cavity // 3 - As above, imaginary response variable doLog	// If true, results are returned in logarithm form, base 10. Default is 0. if (ParamIsDefault(mode)) mode = 0; endif if (ParamIsDefault(doLog)) doLog = 0; endif // Coefficient parameters variable f0 = c[0]			// Natural resonance frequency variable fS = c[1]			// Shifted resonance frequency variable gM = c[2]			// Natural damping rate variable k = c[3]			// Linewidth variable d = c[4]			// Detuning variable aSN = c[5]		// Shot noise power variable aTMax = c[6]		// Maximum technical noise power variable aTMin = c[7]		// Technical noise power offset variable qT = c[8]			// Technical noise peak phase variable q = c[9]			// Quadrature angle variable eDet = c[10]		// Detection efficiency

// Derived parameters variable ksR = fS^2 - f0^2										// Real spring constant variable gOpt = 2*k*(f0^2 - f^2)/(k^2 + d^2 - f^2)					// Optomechanical damping rate variable gTot = gM + gOpt										// Total damping rate variable denomM = (fS^2 - f^2)^2 + f^2*gOpt^2					// Gain denominator magnitude // Intracavity stuff variable GAmR = ((f0^2 - f^2)*(fS^2 - f^2) + f^2*gOpt^2) / denomM	// Real part of the AM gain response to intracavity AM	variable GAmI = (f0^2 - fS^2)*f*gOpt / denomM					// Imaginary part of the AM gain response to intracavity AM	variable GFmR = -(ksR/d)*(k*(fS^2 - f^2) + f^2*GOpt) / denomM	// Real part of the FM gain response to intracavity AM	variable GFmI = -(ksR/d)*(f*(fS^2 - f^2) + f*k*GOpt) / denomM	// Imaginary part of the FM gain response to intracavity AM

// Shot noise stuff if (mode==0 || mode==1) if (mode==1 && mod(q,pi/2)!=0) Abort "Error in omFitMaster: Mode 1 can only be used with AM or FM quadratures (you chose "+num2str(q)+" rad)." endif // Extracavity stuff -- note that the following explicitly ignores the high-frequency filtering of the cavity Lorentzian! variable iRot = 2*k^2/(k^2+d^2-f^2)		// Rotation of drive to same phase quadrature (i.e. AM to AM) variable qRot = 2*k*d/(k^2+d^2-f^2)		// Rotation of drive to out of phase quadrature (i.e. AM to FM) variable MAmSnR	= (GAmR*iRot - 1)*cos(q) 	+ (GFmR*iRot - qRot)*sin(q) 	// Real response to AM shot noise variable MAmSnI	= (GAmI*iRot)*cos(q) 		+ (GFmI*iRot)*sin(q)			// Imaginary response to AM shot noise variable MFmSnR	= (GAmR*qRot)*cos(q) 	+ (iRot + GFmR*qRot - 1)*sin(q)	// Real response to FM shot noise variable MFmSnI 	= (GAmI*qRot)*cos(q) 	+ (GFmI*qRot)*sin(q)			// Imaginary response to FM shot noise variable SSn = MAmSnR^2 + MAmSnI^2 + MFmSnR^2 + MFmSnI^2		// Total power response to shot noise // Technical noise variable DAmT, DQT, SQT if (mode==0) DAmT = aTMin + aTMax*cos(qT)		// Technical noise power in AM DQT = aTMin + aTMax*cos(qT-q)		// Technical noise power at this quadrature SQT = max(0,DQT - DAmT)			// Technical noise contribution from non-AM at this quadrature and frequency // Can not be less than 0 else	// mode 1 wave techAm, techFm DAmT = techAm(f) if (mod(q,pi)==0) SQT = 0 else		// Calculating FM quadrature SQT = techFm(f) endif endif variable SGTech = ( (GAmR*cos(q) + GFmR*sin(q))^2 + (GAmI*cos(q) + GFmI*sin(q))^2 )*DAmT	// Atomic response to technical drive variable STotal = (SSn + SGTech + SQT)*aSN		// Total noise power STotal = eDet*STotal + aSN*(1-eDet)		// Correct for detection efficiency // Return output if (doLog) return log(STotal) else return STotal endif endif

if (mode >= 2 ) // Network analysis stuff variable rotAmFm = (f^2*d^2)/((k^2 + d^2)^2 + f^2*k^2) variable SAm, SFm if (mode==3) SAm = GAmR + 1		// Real AM to real AM			SFm = GFmR + rotAmFm	// Real AM to real FM		elseif (mode==4) SAm = GAmI			// Real AM to imaginary AM			SFm = GFmI				// Real AM to imaginary FM		else abort "Error in omFitMaster. Mode was "+num2str(mode)+", must be between 0 and 3." endif return (SAm*cos(q) + SFm*sin(q))*aTMax endif

end

//--- // function makeOmPhaseMap(c,m) // //	Makes a theoretical map of power response vs frequency and //	quadrature phase. // //	c is a coefficient wave containing: //		c[0] -- Natural frequency //		c[1] -- Shifted frequency //		c[2] -- Natural damping rate (frequency units) //		c[3] -- Cavity half-linewidth (frequency units) //		c[4] -- Cavity detuning (frequency units) //		c[5] -- Shot noise level //		c[6] -- Maximum technical noise level (units of noise*frequency^2) //		c[7] -- Minimum technical noise level (units of noise*frequency^2) //		c[8] -- Noise peak quadrature angle (radians) // 		c[9] -- Detection efficiency //	Here the technical noise is assumed to decay as 1/f^2 in all //	quadratures. // //	m is the destination wave. Frequency is indexed by rows, //	and quadrature phase by columns. The axes must be properly //	scaled before calling makeOmPhaseMap. //--- function makeOmPhaseMap(c,m) wave c, m	// Make the coefficient wave for omFitMaster make /n=11 cTemp cTemp[0,8] = c[p] cTemp[10] = c[9] variable j, i	for (j=0; j<DimSize(m,1); j+=1)	// Loop over quadrature angles cTemp[9] = DimOffset(m,1) + j*DimDelta(m,1)		// This is the quadrature angle m[][j] = omFitMaster(cTemp,x,mode=0) endfor killwaves/z cTemp end

//--- // function omResponse(c,f) : fitfunc // //	Fitting function for linear optomechanical power response //	Coefficients are: //		c[0] -- Natural frequency //		c[1] -- Shifted frequency //		c[2] -- Natural damping rate (frequency units) //		c[3] -- Cavity half-linewidth (frequency units) //		c[4] -- Cavity detuning (frequency units) //		c[5] -- Shot noise level //		c[6] -- Maximum technical noise level //		c[7] -- Minimum technical noise level //		c[8] -- Noise peak quadrature angle (radians) //		c[9] -- Detection quadrature angle (radians) // 		c[10] -- Detection efficiency //--- function omResponse(c,f) : fitfunc wave c	variable f	//CurveFitDialog/ //CurveFitDialog/ Coefficients 11 //CurveFitDialog/ c[0] -- f0	//CurveFitDialog/ c[1] -- fS	//CurveFitDialog/ c[2] -- G	//CurveFitDialog/ c[3] -- k	//CurveFitDialog/ c[4] -- d	//CurveFitDialog/ c[5] -- ASN //CurveFitDialog/ c[6] -- ATH //CurveFitDialog/ c[7] -- ATL //CurveFitDialog/ c[8] -- phiT //CurveFitDialog/ c[9] -- phi //CurveFitDialog/ c[10] -- eps return omFitMaster(c,f,mode=0) end

//--- // function logOmResponse(c,f) : fitfunc // //	Fitting function for linear optomechanical power response, //	returning log of power //	Coefficients are: //		c[0] -- Natural frequency //		c[1] -- Shifted frequency //		c[2] -- Natural damping rate (frequency units) //		c[3] -- Cavity half-linewidth (frequency units) //		c[4] -- Cavity detuning (frequency units) //		c[5] -- Maximum technical noise level (relative to shot noise) //		c[6] -- Minimum technical noise level (relative to shot noise //		c[7] -- Noise peak quadrature angle (radians) //		c[8] -- Detection quadrature angle (radians) // 		c[9] -- Detection efficiency //--- function logOmResponse(c,f) : fitfunc	wave c	variable f	//CurveFitDialog/	//CurveFitDialog/ Coefficients 11	//CurveFitDialog/ c[0] -- f0	//CurveFitDialog/ c[1] -- fS	//CurveFitDialog/ c[2] -- G	//CurveFitDialog/ c[3] -- k	//CurveFitDialog/ c[4] -- d	//CurveFitDialog/ c[5] -- ATH	//CurveFitDialog/ c[6] -- ATL	//CurveFitDialog/ c[7] -- phiT	//CurveFitDialog/ c[8] -- phi	//CurveFitDialog/ c[9] -- eps	// Coefficient wave for omFitMaster as shot noise level = 1	make /free/n=11 c2	c2[0,4] = c[p]	c2[5] = 1	c2[6,10] = c[p-1]	return log(omFitMaster(c2,f,mode=0)) end

//--- //function Fano(w,t) : fitfunc //	Fitting Function for shot-noise driven Fano resonance with //	finite detection efficiency. //--- function fano(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = amp //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = w0	//CurveFitDialog/ w[4] = gamma return w[0]*((1-w[1])+w[1]*((t^2-w[3]^2)^2+t^2*w[4]^2)/((t^2-w[2]^2)^2+t^2*w[4]^2))

end //return to omFunctionList

//--- //function logFano(w,t) : fitfunc //	Fitting Function for Log of Fano // //--- function logFano(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = noise //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = w0	//CurveFitDialog/ w[4] = gamma return Log((1-w[1])+w[1]*((t^2-w[3]^2)^2+t^2*w[4]^2)/((t^2-w[2]^2)^2+t^2*w[4]^2)+w[0])

end //return to omFunctionList

//--- //function BreitWigner(w,t) : fitfunc //	Fitting Function for Breit - Wigner resonance // //--- function BreitWigner(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 4 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0	//CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = gamma return w[1] + w[0]*w[2]^2*w[3]^2/( (t^2 - w[2]^2)^2 + t^2 * w[3]^2) End //--- //function myLor(w,t) : fitfunc //	Fitting Function for Lorenztain with independent amplitude and width // //--- Function myLor(w,t) : FitFunc Wave w	Variable t

//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog. //CurveFitDialog/ Equation: //CurveFitDialog/ f(t) = y0 + Amp* (gamma/2)^2 / ( (t - x0)^2 + (gamma/2)^2 ) //CurveFitDialog/ End of Equation //CurveFitDialog/ Independent Variables 1 //CurveFitDialog/ t	//CurveFitDialog/ Coefficients 4 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0	//CurveFitDialog/ w[2] = x0	//CurveFitDialog/ w[3] = gamma

return w[1] + w[0]* (w[3]/2)^2 / ( (t - w[2])^2 + (w[3]/2)^2 ) End

//function BreitWigner2(w,t) : fitfunc //	Fitting Function for Breit - Wigner resonance with nonzero center frequency // //--- function BreitWigner2(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0	//CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = wOffset return w[1] + w[0]*(w[2]-w[4])^2*w[3]^2/( ((t-w[4])^2 - (w[2]-w[4])^2)^2 + (t-w[4])^2 * w[3]^2) End

//function FourBreitWigner(w,t) : fitfunc //	Fitting Function for Breit - Wigner resonance with nonzero center frequency // //--- function FourBreitWigner(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = Amp //CurveFitDialog/ w[1] = y0	//CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = spacing variable ws0 = w[2], ws1 = w[2]+w[4], ws2 = w[2] + 2*w[4], ws3 = w[2] - w[4] variable bw0 = ws0^2*w[3]^2/( (t^2 - ws0^2)^2 + t^2 * w[3]^2) variable bw1 = ws1^2*w[3]^2/( (t^2 - ws1^2)^2 + t^2 * w[3]^2) variable bw2 = ws2^2*w[3]^2/( (t^2 - ws2^2)^2 + t^2 * w[3]^2) variable bw3 = ws3^2*w[3]^2/( (t^2 - ws3^2)^2 + t^2 * w[3]^2)

return w[1] + w[0]*(bw0 + 2*bw1 + 1.5*bw2 + 0.5*bw3)/5

End //-- //--- //function BothBreitWigner(w,t) : fitfunc //	Fitting Function for Breit - Wigner resonance // //--- function BothBreitWigner(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 6 //CurveFitDialog/ w[0] = AmpRed //CurveFitDialog/ w[1] = y0	//CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = carrier //CurveFitDialog/ w[5] = AmpBlue variable AmpRed = w[0] variable wsRed = w[2] variable carrier = w[4] variable wsBlue = 2*carrier - wsRed variable AmpBlue = w[5] return w[1] + AmpRed*wsRed^2*w[3]^2/( (t^2 - wsRed^2)^2 + t^2 * w[3]^2) + AmpBlue*wsBlue^2*w[3]^2/( (t^2 - wsBlue^2)^2 + t^2 * w[3]^2) End //function BothFourBreitWigner(w,t) : fitfunc //	Fitting Function for Breit - Wigner resonance with nonzero center frequency // //--- function BothFourBreitWigner(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = AmpRed //CurveFitDialog/ w[1] = y0	//CurveFitDialog/ w[2] = ws	//CurveFitDialog/ w[3] = gamma //CurveFitDialog/ w[4] = spacing //CurveFitDialog/ w[5] = carrier //CurveFitDialog/ w[6] = AmpBlue variable AmpRed = w[0] variable wsRed = w[2] variable carrier = w[5] variable wsBlue = 2*carrier - wsRed variable AmpBlue = w[6] variable ws0Red = wsRed, 			ws0Blue = wsBlue variable ws1Red = wsRed+w[4], 		ws1Blue = wsBlue - w[4] variable ws2Red = wsRed + 2*w[4], 	ws2Blue = wsBlue - 2*w[4] variable ws3Red = wsRed - w[4],		ws3Blue = wsBlue + w[4] variable bw0Red = ws0Red^2*w[3]^2/( (t^2 - ws0Red^2)^2 + t^2 * w[3]^2) variable bw1Red = ws1Red^2*w[3]^2/( (t^2 - ws1Red^2)^2 + t^2 * w[3]^2) variable bw2Red = ws2Red^2*w[3]^2/( (t^2 - ws2Red^2)^2 + t^2 * w[3]^2) variable bw3Red = ws3Red^2*w[3]^2/( (t^2 - ws3Red^2)^2 + t^2 * w[3]^2)

variable bw0Blue = ws0Blue^2*w[3]^2/( (t^2 - ws0Blue^2)^2 + t^2 * w[3]^2) variable bw1Blue = ws1Blue^2*w[3]^2/( (t^2 - ws1Blue^2)^2 + t^2 * w[3]^2) variable bw2Blue = ws2Blue^2*w[3]^2/( (t^2 - ws2Blue^2)^2 + t^2 * w[3]^2) variable bw3Blue = ws3Blue^2*w[3]^2/( (t^2 - ws3Blue^2)^2 + t^2 * w[3]^2) variable amp0 = 1, amp1 = 2, amp2 = 1.4, amp3 = 0.6; variable ampSum = amp0 + amp1 + amp2 + amp3

return w[1] + AmpRed*(amp0*bw0Red + amp1*bw1Red + amp2*bw2Red + amp3*bw3Red)/ampSum + AmpBlue*(amp0*bw0Blue + amp1*bw1Blue + amp2*bw2Blue + amp3*bw3Blue)/ampSum

End //- //--- function totalAmResponse5(inx) variable inx wave w = w_coef nvar omTfreq variable t = omTfreq //CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs	//CurveFitDialog/ w[3] = f0	//CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta variable deltanoise = w[7]

variable ReGNum, ImGNum, GDenom, ReG, ImG, A, B, KD, k, d, fs, f0, gM	variable techPwr wave techwaveAM fs = w[2];	f0 = w[3];	k = w[5];		d = w[6]+inx;		gM = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt ReGNum = (fs^2 - f0^2)*(t^2 - fs^2)			// Numerator of real part of gain function ImGNum = t*gTot*(fs^2-f0^2)				// Numerator of imaginary part of gain function GDenom = (t^2 - fs^2)^2 + t^2*gTot^2		// Denominator (real) of gain function ReG = ReGNum/GDenom ImG = ImGNum/GDenom A = 2*k^2/(k^2 + d^2)						// Ratio of 2*kappa^2 to kappa^2 + Delta^2 B = (k^2-d^2)/(k^2+d^2)					// Unperturbed reflected field // Real part of the shot noise response // Corrected by NB Feb 28, 2011 variable ReSn = ( ReG*A + B )^2  +  ( ImG*A*d/k )^2 variable ImSn = ( ImG*A )^2         +   ( (ReG+1)*A*d/k )^2 variable Sam if(w[0]==0) Sam = (ReG+1)^2 + ImG^2 else techPwr = (techWaveAM(t)/w[0] - 1)/w[1] variable Tech = techPwr*( (ReG+1)^2 + ImG^2 )

Sam= w[0]*( w[1]*(ReSn + ImSn + Tech) + 1 - w[1] )		// Total AM optomechanical response endif return 1/sqrt(2*pi)/(DeltaNoise)*Exp(-(inx^2)/(2*DeltaNoise^2))*Sam End //return to omFunctionList

//return to omFunctionList //--- //function totalFmResponse4(w,t) : fitfunc //	Fitting Function for shot-noise driven BW resonance with //	finite detection efficiency and optomechanical damping. //	Corrected for phase error in Gain //--- function totalFmResponse4(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs	//CurveFitDialog/ w[3] = f0	//CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta variable ReGNum, ImGNum, GDenom, ReG, ImG, A, B, KD, k, d, fs, f0, gM	variable techPwr wave techwaveAM wave techwaveFM fs = w[2];	f0 = w[3];	k = w[5];		d = w[6];		gM = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt ReGNum = (fs^2 - f0^2)*(k*(fs^2 - t^2) + t^2*gtot)			// Numerator of real part of gain function ImGNum = t*(fs^2-f0^2)*(k*gtot - (fs^2-t^2)) 				// Numerator of imaginary part of gain function GDenom = d*(fs^2 - t^2)^2 + d*t^2*gTot^2		// Denominator (real) of gain function ReG = ReGNum/GDenom ImG = ImGNum/GDenom A = 2*k^2/(k^2 + d^2)						// Ratio of 2*kappa^2 to kappa^2 + Delta^2

variable ReSn = 2*A*(ReG^2 + ImG^2 + d/k*ReG) variable aToA = ((k^2 + d^2)^2 + t^2*k^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) // AM to AM transduction variable aToF = (t^2*d^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) variable ratioFA = aToF/aToA if(w[0]==0) return ReG^2 + ImG^2 + ratioFA else techPwr = (techWaveAM(t)/w[0] - 1)/w[1] variable Tech = techPwr*(ReG^2 + ImG^2)

return w[0]*( w[1]*(ReSn+1+Tech) + 1 - w[1] ) + techwaveFM(t)	-w[0]	// Total FM optomechanical response endif End //return to omFunctionList function fitavgFmResp5(w,t) : fitfunc wave w	variable t

//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs	//CurveFitDialog/ w[3] = f0	//CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta //CurveFitDialog/ w[7] = DeltaNoise duplicate/o w wtemp wtemp=w variable deltaNoise = w[7] variable /g omTfreq = t	return Integrate1d(totalFmResponse5,-4*deltaNoise,4*deltaNoise,0,100) end

function totalFmResponse5(inx) variable inx wave w = w_coef nvar omTfreq variable t = omTfreq //CurveFitDialog/ //CurveFitDialog/ Coefficients 7 //CurveFitDialog/ w[0] = snLevel //CurveFitDialog/ w[1] = eps //CurveFitDialog/ w[2] = fs	//CurveFitDialog/ w[3] = f0	//CurveFitDialog/ w[4] = gammaM //CurveFitDialog/ w[5] = kappa //CurveFitDialog/ w[6] = Delta variable deltaNoise = w[7] variable ReGNum, ImGNum, GDenom, ReG, ImG, A, B, KD, k, d, fs, f0, gM	variable techPwr wave techwaveAM wave techwaveFM fs = w[2];	f0 = w[3];	k = w[5];		d = w[6] + inx;		gM = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt ReGNum = (fs^2 - f0^2)*(k*(fs^2 - t^2) + t^2*gtot)			// Numerator of real part of gain function ImGNum = t*(fs^2-f0^2)*(k*gtot - (fs^2-t^2)) 				// Numerator of imaginary part of gain function GDenom = d*(fs^2 - t^2)^2 + d*t^2*gTot^2		// Denominator (real) of gain function ReG = ReGNum/GDenom ImG = ImGNum/GDenom A = 2*k^2/(k^2 + d^2)						// Ratio of 2*kappa^2 to kappa^2 + Delta^2 variable aToA = ((k^2 + d^2)^2 + t^2*k^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) // AM to AM transduction variable aToF = (t^2*d^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) variable ratioFA = aToF/aToA

variable ReSn = 2*A*(ReG^2 + ImG^2 + d/k*ReG) variable Sfm if(w[0]==0) Sfm = ReG^2 + ImG^2 + ratioFA else techPwr = (techWaveAM(t)/w[0] - 1)/w[1] variable Tech = techPwr*(ReG^2 + ImG^2)

Sfm = w[0]*( w[1]*(ReSn+1+Tech) + 1 - w[1] ) + techwaveFM(t)	-w[0]	// Total FM optomechanical response endif

return 1/sqrt(2*pi)/(DeltaNoise)*Exp(-(inx^2)/(2*DeltaNoise^2))*Sfm End //return to omFunctionList //--- //function AMresponseAMdrive(w,t) : fitfunc //	Fitting Function for AM technical noise driven data //	Plots the response in the AM quadrature of //	(Amplitude_Atoms - Amplitude_SN)/(Amplitude_Empty - Amplitude_SN) //--- function AMresponseAMdrive(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt Variable ResponseNum = ((fs^2-t^2)*(f0^2-t^2)+(t*gTot)^2)^2+(t*gTot*(fs^2-f0^2))^2 Variable ResponseDenom = ((fs^2-t^2)^2+(t*gTot)^2)^2 Return ResponseNum/ResponseDenom

End //return to omFunctionList

//--- //function FMresponseAMdrive(w,t) : fitfunc //	Fitting Function for AM technical noise driven data //	Plots the response in the FM quadrature of //	(Amplitude_Atoms - Amplitude_SN)/(Amplitude_Empty - Amplitude_SN) //--- function FMresponseAMdrive(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt Variable RealNum = d^2*t*((fs^2-t^2)+(t*gTot)^2)+k*(d^2+k^2)*t*gTot*(fs^2-f0^2) Variable ImagNum = -k*(d^2+k^2)*(fs^2-f0^2)*(fs^2-t^2) Variable ResponseDenom = d*(d^2+k^2)*((fs^2-t^2)^2+(t*gTot)^2) variable aToA = ((k^2 + d^2)^2 + t^2*k^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) // AM to AM transduction variable aToF = (t^2*d^2)/((k^2+d^2-t^2)^2 + 4*t^2*k^2) variable ratioFA = aToF/aToA Return (RealNum^2+ImagNum^2)/(ResponseDenom)^2 + ratioFA

End //return to omFunctionList

//--- //function AMIPresponseAMdrive(w,t) : fitfunc //	Fitting Function for the normalized in-phase response in the //	AM quadrature to an AM drive //	The function fits to the following curve //	(Note: Amp = Amplitude and SN = Shot-noise) //	(Amp_Atoms-Amp_SN)/(Amp_EmptyCavity-Amp_SN) //--- function AMIPresponseAMdrive(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM+gOpt Variable ReG = ((f0^2-fs^2)*(fs^2-t^2))/((fs^2-t^2)^2+(t*gTot)^2) Return 1+ReG

End //return to omFunctionList

//--- //function AMOPresponseAMdrive(w,t) : fitfunc //	Fitting Function for the normalized out-of-phase response in the //	AM quadrature to an AM drive //	The function fits to the following curve //	(Note: Amp = Amplitude and SN = Shot-noise) //	Log((Amp_Atoms-Amp_SN)/(Amp_EmptyCavity-Amp_SN)) //--- function AMOPresponseAMdrive(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt Variable ImG = (f0^2-fs^2)*(t*gTot)/((fs^2-t^2)^2+(t*gTot)^2) Return -ImG

End //return to omFunctionList

//--- //function FMresponseAMdrive(w,t) : fitfunc //	Fitting Function for the normalized in-phase response in the //	FM quadrature to an AM drive //	The function fits to the following curve //	(Note: Amp = Amplitude and SN = Shot-noise) //	(Amp_Atoms-Amp_SN)/(Amp_EmptyCavityAM-Amp_SN) //--- function FMIpResponseAmDrive(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt

Variable ReG = ((f0^2-fs^2)*(fs^2-t^2))/((fs^2-t^2)^2+(t*gTot)^2) variable ratioReFA = d*t/sqrt((k^2+d^2)^2 + t^2*k^2) // May need a - sign prefactor depending on Fourier convention Return -k*ReG/d + ratioReFA

End //return to omFunctionList

//--- //function FMresponseAMdrive(w,t) : fitfunc //	Fitting Function for the normalized out-of-phase response in the //	FM quadrature to an AM drive //	The function fits to the following curve //	(Note: Amp = Amplitude and SN = Shot-noise) //	Log((Amp_Atoms-Amp_SN)/(Amp_EmptyCavityAM-Amp_SN)) //--- function FmOpResponseAMdrive(w,t) wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 8 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta Variable fs = w[0]; Variable f0 = w[1]; Variable gM = w[2]; Variable k = w[3]; Variable d = w[4];

variable gOpt = 2*k*(f0^2-t^2)/(k^2 + d^2 - t^2) variable gTot = gM + gOpt // This is what Thierry calculates Variable ImG = (f0^2-fs^2)*(t*gTot)/((fs^2-t^2)^2+(t*gTot)^2) // This is what I think it may be, but is from what I used to calculate the SN response //	this has extra term from -i*w in k - i*w in numerator: //	Variable ImG = t*(fs^2 - f0^2)*(gTot - (fs^2 - t^2)/k) / ( (fs^2- t^2)^2 + t^2*gTot^2) Return k*ImG/d

End //return to omFunctionList

function AllphaseAllquadResponseAmDrive(w,t) wave w	variable t	nvar npoints //CurveFitDialog/ //CurveFitDialog/ Coefficients 5 //CurveFitDialog/ w[0] = fs	//CurveFitDialog/ w[1] = f0	//CurveFitDialog/ w[2] = gammaM //CurveFitDialog/ w[3] = kappa //CurveFitDialog/ w[4] = delta wave rAmIp variable deltaFreq = npoints*deltax(rAmIp) if(t< (leftx(rAmIp) + deltaFreq)) return AMIPresponseAMdrive(w,t) elseif(t< (leftx(rAmIp) + 2*deltaFreq)) return AMOPresponseAMdrive(w,t-deltaFreq) elseif(t< (leftx(rAmIp) + 3*deltaFreq)) return FMIPresponseAMdrive(w,t-2*deltaFreq) else return FMOPresponseAMdrive(w,t-3*deltaFreq) endif end

function makeRelief34(zmin,zmax) variable zmin, zmax

colortab2wave relief

duplicate/o m_colors ctab setscale /i x,zmin,zmax,ctab

variable i for (i=126; i>=0; i-=2) deletepoints (i),1,ctab endfor

end

//return to omFunctionList

//--- //function SSresponse (w,t) : fitfunc //	Fitting Function for the single-sideband on-resonance om //	response // //--- function SSresponse(w,t) : fitfunc wave w	variable t	//CurveFitDialog/ //CurveFitDialog/ Coefficients 10 //CurveFitDialog/ w[0] = gom //CurveFitDialog/ w[1] = wm	//CurveFitDialog/ w[2] = Gamma_m //CurveFitDialog/ w[3] = p_th //CurveFitDialog/ w[4] = SN_Amplitude //CurveFitDialog/ w[5] = Offset //CurveFitDialog/ w[6] = Het_Det_Efficiency //CurveFitDialog/ w[7] = Kappa //CurveFitDialog/ w[8] = w_cavity //CurveFitDialog/ w[9] = nbar variable gom = w[0] variable wm = w[1] variable Gm = w[2] variable pth = w[3] variable SNamp = w[4]	// Ideally, this should equal P_LO variable SNoffset = w[5]	// Ideally, this should equal 1 variable EffHet = w[6] variable kappa = w[7] variable fc = w[8]			// frequency of cavity resonance variable nbar = w[9] variable Copt = 4*nbar*gom^2 / kappa / Gm	variable denom, numerator denom = (wm^2+(Gm/2)^2-(t-fc)^2)^2 + ((t-fc)*Gm)^2 numerator = 2 * wm * (Copt*wm-(t-fc)) + (wm^2 + (Gm/2)^2+(t-fc)^2) * (2*pth+1) return SNamp *(SNoffset + EffHet*Copt*Gm^2 * numerator / denom) End

//return to omFunctionList