(*The basic function, derived from Voinova et al. Phys. Scr. 1999, 59, 391-396
β is the complex response function (see below)
μfilm is the film's elastic shear modulus
ηfilm is the film's dynamic (shear) viscosity
ρfilm is the film's density
hfilm is the film's thickness
ηbulk is the fluid's dynamic viscosity
ρbulk is the fluid's density
I is the square root of -1
Harm is the harmonic number

QCM-D Frequency reponse: Δf = Im[β/(2π ρq hq)] where ρq is the density of the quartz resonator and hq is its thickness
QCM-D Dissipation response: Δd = -Re[β/(π Harm f1 ρq hq)] where f1 is the fundamental resonant frequency of the quartz resonator
*)
$Assumptions=μfilm∈Reals&&ηfilm∈Reals&&ρfilm∈Reals&&hfilm∈Reals&&ηbulk∈Reals&&ρbulk∈Reals&&f1∈Reals&&hq∈Reals&&ρq∈Reals&&Harm∈Reals&&μfilm>0&&ηfilm>0&&ρfilm>0&&hfilm>0&&ηbulk>0&&ρbulk>0&&f1>0&&hq>0&&ρq>0&&Harm>0&&Harm∈Integers
β[μfilm_,ηfilm_,ρfilm_,hfilm_,ηbulk_,ρbulk_]= FullSimplify[Simplify[Refine[κ1 ξ1 (1-A Exp[2 ξ1 hfilm])/(1+A Exp[2 ξ1 hfilm])//.{κ1->ηfilm-I μfilm/(2π Harm f1),κf->ηbulk,ξ1->Sqrt[-ρfilm (2π Harm f1)^2/(μfilm+I ηfilm (2π Harm f1))],ξf->Sqrt[I ρbulk (2 π Harm f1)/ηbulk],A->(κ1 ξ1+κf ξf)/(κ1 ξ1-κf ξf)}]]]

(*Derivative of the beta function with respect to film density ∂β/∂ρfilm.
This needs to be done in stages for Mathematica to be able to process it in a sensible length of time
*)
$Assumptions = True
FullSimplfy[Simplify[Refine[D[β[μfilm,ηfilm,ρfilm,hfilm,ηbulk,ρbulk],ρfilm],Assumptions->μfilm∈Reals&&ηfilm∈Reals&&ρfilm∈Reals&&hfilm∈Reals&&ηbulk∈Reals&&ρbulk∈Reals&&f1∈Reals&&hq∈Reals&&ρq∈Reals&&Harm∈Reals&&μfilm>0&&ηfilm>0&&ρfilm>0&&hfilm>0&&ηbulk>0&&ρbulk>0&&f1>0&&hq>0&&ρq>0&&Harm>0&&Harm∈Integers]]](*The assumptions here are the same as those included below using the $Assumptions tag*)

(*Below is included for comparison and can be deleted later*)

(-I f1 Harm (Sqrt[2 π]) (Sqrt[-2 I f1 Harm π ηfilm-μfilm]) (Sqrt[ρfilm]) ρbulk Cosh[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[-2 I f1 Harm π ηfilm-μfilm]))]+(-2 f1 Harm π ηfilm+I μfilm) ρfilm (Sqrt[((I f1 Harm ρbulk)/(ηbulk))]) Sinh[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[-2 I f1 Harm π ηfilm-μfilm]))])/((Sqrt[-2 I f1 Harm π ηfilm-μfilm]) (Sqrt[ρfilm]) (Sqrt[((I f1 Harm ρbulk)/(ηbulk))]) Cosh[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[-2 I f1 Harm π ηfilm-μfilm]))]+f1 Harm (Sqrt[2 π]) ρbulk Sinh[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[-2 I f1 Harm π ηfilm-μfilm]))])

(*Assumptions are that all material parameters are positive real numbers and that the harmonic number is an integer*)$Assumptions=μfilm∈Reals&&ηfilm∈Reals&&ρfilm∈Reals&&hfilm∈Reals&&ηbulk∈Reals&&ρbulk∈Reals&&f1∈Reals&&hq∈Reals&&ρq∈Reals&&Harm∈Reals&&μfilm>0&&ηfilm>0&&ρfilm>0&&hfilm>0&&ηbulk>0&&ρbulk>0&&f1>0&&hq>0&&ρq>0&&Harm>0&&Harm∈Integers

(*Derivative of the beta function with respect to film shear modulus ∂β/∂μfilm*)
FullSimplify[D[β[μfilm,ηfilm,ρfilm,hfilm,ηbulk,ρbulk],μfilm]]

(f1 Harm (4 I (f1^2) hfilm (Harm^2) π μfilm ((ρfilm)^(7/2))+8 (f1^3) hfilm (Harm^3) (π^2) ((ρfilm)^(5/2)) (-ηfilm ρfilm+ηbulk ρbulk)+(2+2 I) (Sqrt[π]) (μfilm (Sqrt[(f1^3) (Harm^3) ηbulk ((ρfilm)^5) ρbulk])+2 I π ηfilm (Sqrt[(f1^5) (Harm^5) ηbulk ((ρfilm)^5) ρbulk]))-(2+2 I) (Sqrt[π]) (μfilm (Sqrt[(f1^3) (Harm^3) ηbulk ((ρfilm)^5) ρbulk])+2 I π ηfilm (Sqrt[(f1^5) (Harm^5) ηbulk ((ρfilm)^5) ρbulk])) Cos[((4 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]+f1 Harm (Sqrt[2 I f1 Harm π ηfilm+μfilm]) ((ρfilm)^2) (-I μfilm ρfilm+2 f1 Harm π (ηfilm ρfilm+ηbulk ρbulk)) Sin[((4 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]))/(4 (((2 f1 Harm π ηfilm-I μfilm))^2) (Sqrt[ρfilm]) (((I f1 Harm ρfilm Cos[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]+(1+I) (Sqrt[π]) (Sqrt[-(((f1^3) (Harm^3) ηbulk ρfilm ρbulk)/(2 I f1 Harm π ηfilm+μfilm))]) Sinh[2 f1 hfilm Harm π (Sqrt[-((ρfilm)/(2 I f1 Harm π ηfilm+μfilm))])]))^2))

(*Derivative of the beta function with respect to film viscosity ∂β/∂ηfilm*)
FullSimplify[D[β[μfilm,ηfilm,ρfilm,hfilm,ηbulk,ρbulk],ηfilm]]

((f1^2) (Harm^2) π ((ρfilm)^2) (-4 (f1^2) hfilm (Harm^2) π μfilm (Sqrt[2 I f1 Harm π ηfilm+μfilm]) ((ρfilm)^(3/2)) ρbulk+8 I (f1^3) hfilm (Harm^3) (π^2) ρbulk (-ηfilm (Sqrt[2 I f1 Harm π ηfilm+μfilm]) ((ρfilm)^(3/2))+ηbulk (Sqrt[(2 I f1 Harm π ηfilm+μfilm) ρfilm]) ρbulk)-(2-2 I) (Sqrt[π]) (μfilm (Sqrt[(f1^3) (Harm^3) ηbulk (2 I f1 Harm π ηfilm+μfilm) ρfilm ((ρbulk)^3)])+2 I π ηfilm (Sqrt[(f1^5) (Harm^5) ηbulk (2 I f1 Harm π ηfilm+μfilm) ρfilm ((ρbulk)^3)]))+(2+2 I) (Sqrt[π]) (-I μfilm (Sqrt[(f1^3) (Harm^3) ηbulk (2 I f1 Harm π ηfilm+μfilm) ρfilm ((ρbulk)^3)])+2 π ηfilm (Sqrt[(f1^5) (Harm^5) ηbulk (2 I f1 Harm π ηfilm+μfilm) ρfilm ((ρbulk)^3)])) Cos[((4 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]-f1 Harm (2 f1 Harm π ηfilm-I μfilm) ρbulk (-I μfilm ρfilm+2 f1 Harm π (ηfilm ρfilm+ηbulk ρbulk)) Sin[((4 f1 hfilm Harm π)/(Sqrt[((2 I f1 Harm π ηfilm+μfilm)/(ρfilm))]))]))/(2 (((2 f1 Harm π ηfilm-I μfilm))^2) (Sqrt[(2 I f1 Harm π ηfilm+μfilm) ρfilm]) ρbulk (((I f1 Harm ρfilm Cos[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]+(1+I) (Sqrt[π]) (Sqrt[-(((f1^3) (Harm^3) ηbulk ρfilm ρbulk)/(2 I f1 Harm π ηfilm+μfilm))]) Sinh[2 f1 hfilm Harm π (Sqrt[-((ρfilm)/(2 I f1 Harm π ηfilm+μfilm))])]))^2))

(*Derivative of the beta function with respect to film thickness ∂β/∂hfilm*)
FullSimplify[D[β[μfilm,ηfilm,ρfilm,hfilm,ηbulk,ρbulk],hfilm]]

(2 (f1^3) (Harm^3) π ρfilm (-I μfilm ρfilm+2 f1 Harm π (ηfilm ρfilm-ηbulk ρbulk)))/(((((I f1 Harm ρfilm Cos[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))])/(Sqrt[-((ρfilm)/(2 I f1 Harm π ηfilm+μfilm))]))+(1+I) (Sqrt[π]) (Sqrt[(f1^3) (Harm^3) ηbulk ρbulk]) Sinh[2 f1 hfilm Harm π (Sqrt[-((ρfilm)/(2 I f1 Harm π ηfilm+μfilm))])]))^2)

(*Derivative of the beta function with respect to fluid viscosity ∂β/∂ηbulk*)
FullSimplify[D[β[μfilm,ηfilm,ρfilm,hfilm,ηbulk,ρbulk],ηbulk]]

-((((1/2)+(I/2)) (Sqrt[π]) (Sqrt[((f1 Harm ρbulk)/(ηbulk))]))/(((Cos[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]-(((1-I) f1 Harm (Sqrt[π]) ηbulk ρbulk Sin[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))])/(Sqrt[f1 Harm ηbulk (2 I f1 Harm π ηfilm+μfilm) ρfilm ρbulk]))))^2))

(*Derivative of the beta function with respect to fluid density ∂β/∂ρbulk*)
FullSimplify[D[β[μfilm,ηfilm,ρfilm,hfilm,ηbulk,ρbulk],ρbulk]]

-((I f1 Harm (Sqrt[(π/2)]))/((Sqrt[((I f1 Harm ρbulk)/(ηbulk))]) (((Cos[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))]-(((1-I) f1 Harm (Sqrt[π]) ηbulk ρbulk Sin[((2 f1 hfilm Harm π (Sqrt[ρfilm]))/(Sqrt[2 I f1 Harm π ηfilm+μfilm]))])/(Sqrt[f1 Harm ηbulk (2 I f1 Harm π ηfilm+μfilm) ρfilm ρbulk]))))^2)))

(*Christopher Blanford August 2015*)