(*^ ::[ frontEndVersion = "Microsoft Windows Mathematica Notebook Front End Version 2.2"; microsoftWindowsStandardFontEncoding; fontset = title, "Times", 24, L1, center, nohscroll, bold; fontset = subtitle, "Times", 18, L1, center, nohscroll, bold; fontset = subsubtitle, "Times", 14, L1, center, nohscroll, italic; fontset = section, "Times", 18, L1, nohscroll, bold, grayBox; fontset = subsection, "Times", 14, L1, nohscroll, bold, blackBox; fontset = subsubsection, "Times", 12, L1, nohscroll, bold, whiteBox; fontset = text, "Times New Roman", 12, L1, nohscroll; fontset = smalltext, "Times", 10, L1, nohscroll; fontset = input, "Courier", 12, L1, nowordwrap; fontset = output, "Courier", 12, L-5, nowordwrap; fontset = message, "Courier", 12, L1, nowordwrap; fontset = print, "Courier", 12, L1, nowordwrap; fontset = info, "Courier", 12, L1, nowordwrap; fontset = postscript, "Courier", 12, L1, nowordwrap; fontset = name, "Times", 10, L1, nohscroll, italic; fontset = header, "Times New Roman", 12, L1, nohscroll; fontset = footer, "Times New Roman", 12, L1, center, nohscroll; fontset = help, "Times", 10, L1, nohscroll; fontset = clipboard, "Times New Roman", 12, L1, nohscroll; fontset = completions, "Courier", 12, L1, nohscroll; fontset = graphics, "Courier New", 10, L0, nowordwrap, nohscroll; fontset = special1, "Times New Roman", 12, L1, nohscroll; fontset = special2, "Times New Roman", 12, L1, nohscroll; fontset = special3, "Times New Roman", 12, L1, nohscroll; fontset = special4, "Times New Roman", 12, L1, nohscroll; fontset = special5, "Times New Roman", 12, L1, nohscroll; fontset = leftheader, "Times New Roman", 12, L0, nohscroll, cellOutline; fontset = leftfooter, "Times New Roman", 12, L0, blackBox, cellOutline; fontset = reserved1, "Courier New", 10, L0, nowordwrap, nohscroll;] :[font = subtitle; inactive; preserveAspect; nohscroll; center; ] Mode-Coupling LG Solution for Cochlear Mechanics :[font = subsubtitle; inactive; preserveAspect; nohscroll; center; ] Lloyd Watts CALTECH 116-81 Pasadena CA 91125 lloyd@hobiecat.cs.caltech.edu :[font = text; inactive; preserveAspect; nohscroll; ] This notebook contains commands to implement the Mode-Coupling Liouville-Green solution for cochlear mechanics, as described in the Ph.D. Thesis by Lloyd Watts. :[font = input; preserveAspect; startGroup; nowordwrap; ] (* Generic Commands, load once only *) :[font = input; initialization; preserveAspect; startGroup; nowordwrap; ] *) Off[General::spell]; Off[FindRoot::frmp]; <Identity; display = DisplayFunction->$DisplayFunction; $DefaultFont = {"Helvetica",10}; (* $DefaultFont = {"Times-Roman",10}; *) beeptable = Table[N[Sin[Pi t/0.2]^2 Sin[1500 t],5],{t,0,0.2,.001 0.2}]; beep := ListPlay[beeptable]; (* :[font = message; inactive; preserveAspect; nowordwrap; ] CellArray::obsfn: CellArray is an obsolete function, superseded by {Raster, RasterArray}. :[font = input; initialization; preserveAspect; nowordwrap; ] *) Clear[phase,expintlist,expdbydxlist,continuous,expinterpolation3,db]; phase[y_] := Module[{sum}, Table[ If[i==1,sum=Arg[y[[1]]],sum+=Arg[y[[i]]/y[[i-1]] ] ], {i,Length[y]}] ]; expintlist[x_, y_] := Module[{sum}, Table[ If[i==1,sum=0, (*else*)sum+=(y[[i]]-y[[i-1]])*(x[[i]]-x[[i-1]])/ Log[E,y[[i]]/y[[i-1]]]//N], {i,Length[y]}] ]; expdbydxlist[x_, y_] := Table[ Which[i==1,y[[i]]/(x[[i+1]]-x[[i]])*Log[E,y[[i+1]]/y[[i]]]//N, i>1&&ix[[n]],n++]; lown = Max[n-1,1]; highn = Max[n,2]; If[lown==0 || y[[highn]]==0,0, Exp[Log[E,y[[lown]]] + (clipx-x[[lown]])/(x[[highn]]-x[[lown]])* Log[E,y[[highn]]/y[[lown]]] ] //N] ]; db[y_] := 20 Log[10,Abs[y] ]; (* :[font = input; initialization; preserveAspect; endGroup; endGroup; nowordwrap; ] *) listplot[list_, options___] := ListPlot[Transpose[{xtable,list}],options,Frame->True,PlotJoined->True, AspectRatio->0.8,Axes->False]; dbydxlistalt[list_] := Module[{diff},Table[ Which[i==1,diff=(list[[2]]-list[[1]])//N, i>1&&i1&&iTrue,PlotJoined->True, PlotRange->{{1,pointdens*h+1},All},AspectRatio->0.8,Axes->False, FrameTicks->{Table[{i,(i-1)*h/(pointdens)//N},{i,1,pointdens h+1, (pointdens h )/ 5}], Automatic} ]; verticalplot2[grid_, options___] := MultipleListPlot[Take[Transpose[grid], {takexmin*pointdens+1,takexmax*pointdens+1}],options, Frame->True,PlotJoined->True, PlotRange->{{1,pointdensy*h+1},All},AspectRatio->0.8,Axes->False, FrameTicks->{Table[{i,(i-1)*h/(pointdensy)//N},{i,1,pointdensy h+1, (pointdensy h )/ 5}], Automatic} ]; horizontalplot[grid_, options___] := MultipleListPlot[Transpose[Take[Transpose[grid], {takexmin*pointdens+1,takexmax*pointdens+1}]],options, Frame->True,PlotJoined->True, PlotRange->{All,All},AspectRatio->0.8,Axes->False, FrameTicks->{Table[{pointdens*(i-takexmin)+1,i},{i,takexmin,takexmax}], Automatic} ]; surfaceplot[grid_, options___] := ListPlot3D[Transpose[Take[Transpose[grid], {takexmin*pointdens+1,takexmax*pointdens+1}]],options, BoxRatios->{takexmax-takexmin,h,h},PlotRange->All, Ticks->{Table[{pointdens*(i-takexmin)+1,i},{i,takexmin,takexmax}], Table[{i,(i-1)*h/(pointdens-1)//N}, {i,1,pointdens h,(pointdens-1) h / 5}], Automatic},ClipFill->None ]; laplacecheckplot[grid_, options___] := ListDensityPlot[-Table[If[j > 1 && j < Length[grid[[1]]] && i > 1 && i < Length[grid], Abs[1 + (grid[[i+1,j]] + grid[[i-1,j]] - 2 grid[[i,j]])/ (grid[[i,j+1]] + grid[[i,j-1]] - 2 grid[[i,j]])], 0],{i,1,Length[grid]},{j,1,Length[grid[[1]]]}] ,AspectRatio->1/10,Mesh->False,PlotRange->{-1,0}, FrameTicks->{Table[{pointdens*i,i},{i,0,xmax}], Table[{i,i*h/(pointdens+1)//N},{i,0,(pointdens+1) h,(pointdens+1)h/5}]}, options]; laplacecheck2[grid_] := Table[If[j > 1 && j < Length[grid[[1]]] && i > 1 && i < Length[grid], grid[[i+1,j]] + grid[[i-1,j]] + grid[[i,j+1]] + grid[[i,j-1]] - 4 grid[[i,j]], 0],{i,1,Length[grid]},{j,1,Length[grid[[1]]]}]; gridplot[grid_, options___] := ListDensityPlot[grid,options,AspectRatio->1/10,Mesh->False,PlotRange->{-1,0}, FrameTicks->{Table[{pointdens*i,i},{i,0,xmax}], Table[{i,i*h/(pointdens+1)//N},{i,0,(pointdens+1) h, (pointdens+1) h/5}]}]; sliceplot[list_, options___] := ListPlot[Take[list, {takexmin*pointdens+1,takexmax*pointdens+1}],options, Frame->True,PlotJoined->True, PlotRange->{All,All},AspectRatio->0.8,Axes->False, FrameTicks->{Table[{pointdens*(i-takexmin)+1,i},{i,takexmin,takexmax}], Automatic} ]; verticalintgrid[grid_] := Module[{sum},Table[CompoundExpression[ Do[ If[i==1,sum=0,sum+=(grid[[i,j]]+grid[[i-1,j]])/2//N], {i,1,Length[grid]}]; sum] ,{j,1,Length[grid[[1]]]}]]; argplot[grid_, options___] := ListPlot[Transpose[{Re[grid],Im[grid]}],options,PlotJoined->True, AspectRatio->Automatic,PlotRange->All,Frame->True,Axes->False]; argpointplot[grid_, options___] := Show[ ListPlot[Transpose[{Re[grid],Im[grid]}],PlotJoined->True,PlotRange->All, nodisplay], ListPlot[Transpose[{Re[grid],Im[grid]}],PlotRange->All,nodisplay, PlotStyle->PointSize[0.015]],options,display, AspectRatio->Automatic,PlotRange->All,Frame->True,Axes->False]; listpointplot[list_, options___] := Show[ ListPlot[list,PlotJoined->True,PlotRange->All, nodisplay], ListPlot[list,PlotRange->All,nodisplay, PlotStyle->PointSize[0.015]],options,display, AspectRatio->Automatic,PlotRange->All,Frame->True,Axes->False]; complexlistplot[list_, options___] := Show[ListPlot[Re[list],PlotJoined->True,nodisplay], ListPlot[Im[list],PlotJoined->True,PlotStyle->gray3,nodisplay], options,display,AspectRatio->0.8,Frame->True,Axes->False,PlotRange->All]; complexlistpointplot[list_, options___] := Show[ListPlot[Re[list],PlotJoined->True,nodisplay], ListPlot[Re[list],nodisplay,PlotStyle->PointSize[0.015]], ListPlot[Im[list],PlotJoined->True,PlotStyle->gray3,nodisplay], ListPlot[Im[list],PlotStyle->{PointSize[0.015],gray3},nodisplay], options,display,AspectRatio->0.8,Frame->True,Axes->False,PlotRange->All]; ydb[array_] := Transpose[{Transpose[array][[1]],db[Transpose[array][[2]]]}]; yphase[array_] := Transpose[{Transpose[array][[1]],phase[Transpose[array][[2]]]}]; yre[array_] := Transpose[{Transpose[array][[1]],Re[Transpose[array][[2]]]}]; yim[array_] := Transpose[{Transpose[array][[1]],Im[Transpose[array][[2]]]}]; (******************************************************************** expinterpolation : this program returns a value at a given position, interpolated at a given position from a given sampled function. The sampled function must be given in terms of its magnitude (dB) and phase. the phase must be continuous, i.e. have been unwrapped using a program like phaseclean(). this variation assumes x values are equally spaced. ********************************************************************) expinterpolation[xtable_, maglist_, phaselist_, x_] := Module[{numpoints, dx, posn, magdb1, magdb2, magdb,phase1,phase2, phaseval,returnval}, CompoundExpression[ numpoints = Length[xtable]; dx = (xtable[[numpoints]] - xtable[[1]])/(numpoints-1); posn = (x - xtable[[1]])/dx + 1; If[posn >= numpoints,magdb = maglist[[Floor[posn] ]], CompoundExpression[ magdb1 = maglist[[Floor[posn] ]]; magdb2 = maglist[[Floor[posn]+1 ]]; magdb = magdb1 + (magdb2 - magdb1)*(posn-Floor[posn])]]; If[posn >= numpoints,phaseval = phaselist[[Floor[posn] ]], CompoundExpression[ phase1 = phaselist[[Floor[posn] ]]; phase2 = phaselist[[Floor[posn]+1 ]] ; phaseval = phase1 + (phase2 - phase1)*(posn-Floor[posn])]]; returnval = 10^(magdb/20)*Exp[I phaseval]//N; (* Print["x = ",x," value = ",returnval]; *) Return[returnval] ]]; (******************************************************************** expinterpolation2 : this program returns a value at a given position, interpolated at a given position from a given sampled function. The sampled function must be given in terms of its magnitude (dB) and phase. the phase must be continuous, i.e. have been unwrapped using a program like phaseclean(). this particular variation does not assume equally spaced x values. ********************************************************************) expinterpolation2[xtable_, maglist_, phaselist_, x_] := Module[{magdb1, magdb2, magdb,phase1,phase2, phaseval,returnval, lown, highn, n, clipx}, CompoundExpression[ clipx = Min[x,xtable[[Length[xtable]]] ]; (* check for out-of-range *) n = 1; While[clipx>xtable[[n]],n++]; lown = Max[n-1,1]; highn = Max[n,2]; magdb1 = maglist[[lown]]; magdb2 = maglist[[highn]]; magdb = magdb1 + (magdb2 - magdb1)*(clipx-xtable[[lown]])/ (xtable[[highn]]-xtable[[lown]]); phase1 = phaselist[[lown]]; phase2 = phaselist[[highn]] ; phaseval = phase1 + (phase2 - phase1)*(clipx-xtable[[lown]])/ (xtable[[highn]]-xtable[[lown]]); returnval = 10^(magdb/20)*Exp[I phaseval]//N; (* Print["x = ",x," value = ",returnval]; *) Return[returnval] ]]; (******************************************************************** lininterpolation : this program returns a value at a given position, interpolated at a given position from a given sampled function. does not assume equally spaced x values. ********************************************************************) lininterpolation[xtable_, list_, x_] := Module[{returnval, lown, highn, n, clipx}, CompoundExpression[ clipx = Min[x,xtable[[Length[xtable]]] ]; (* check for out-of-range *) n = 1; While[clipx>xtable[[n]],n++]; lown = Max[n-1,1]; highn = Max[n,2]; returnval = list[[lown]] + (list[[highn]]-list[[lown]])* (clipx-xtable[[lown]])/(xtable[[highn]]-xtable[[lown]])//N; (* Print["x = ",x," value = ",returnval]; *) Return[returnval] ]]; (* may be able to use these somehow to save smaller files *) round[j_] := N[If[Re[j]==0,0, Round[10^6 MantissaExponent[N[Re[j]]][[1]]]/10^6 * 10^MantissaExponent[ N[Re[j]] ][[2]] ] + If[Im[j]==0,0, I*Round[10^6 MantissaExponent[N[Im[j]]][[1]]]/10^6 * 10^MantissaExponent[N[Im[j]]][[2]] ] ]; roundlist[list_] := Table[round[list[[i]]],{i,1,Length[list]}]; roundarray[array_] := Table[roundlist[array[[i]]],{i,1,Length[array]}]; (* :[font = input; preserveAspect; startGroup; nowordwrap; ] (* Physical parameters *) :[font = input; initialization; preserveAspect; endGroup; nowordwrap; ] *) (* PHYSICAL PARAMETERS *) rho = .001; (* density of water in g/mm3 *) xmax = 20; (* length of cochlea in mm *) h = 1; (* height of scalae in mm *) d = 5; (* characteristic length of cochlea in mm *) s = 10 10^6 Exp [-x/d]; (* membrane stiffness/area in g/(s2mm2) *) beta = 2; (* membrane damping/area in g/smm2 *) m = 1.5 10^(-3); (* membrane mass/area in g/mm2 *) fo = 1600 Sqrt[2]//N; (* frequency of input in Hz *) (* COMPUTATIONAL PARAMETERS *) dx = 1/7; (* point spacing in mm for finite-diff *) guess1 = .01-.01 I; (* starting point for root near 0 for LG*) guess2 = .01-(Pi-.1) I//N; (* starting point for root near -I Pi for LG*) maxdeltak = .2; (* maximum step in k for LG *) maxdeltax = 1; (* maximum step in x for LG *) (* DERIVED QUANTITIES *) xtable = Table[x,{x,0,xmax,dx}];(* table of x values *) pointdens = 1/dx; (* :[font = input; preserveAspect; startGroup; Cclosed; nowordwrap; ] (* DEFINE: Finite Difference Method method *) :[font = input; initialization; preserveAspect; endGroup; nowordwrap; ] *) fd[rho_, xmax_, h_, d_, s_, (* inputs, named as above *) beta_, m_, fo_, dx_, (* inputs, named as above *) disp_, (* membrane displacement (1-D list) *) press_ (* fluid pressure (2-D list) *) ] := Module[{ (* LOCAL VARIABLES *) a, (* A matrices *) b, (* B matrices *) c, (* C matrices *) p, (* P vectors (differential pressure) *) y, (* Y membrane admittance *) q, (* Q vector *) nx, (* number of grid points in X direction *) ny, (* number of grid points in Y direction *) wo (* angular frequency in radians/s *) }, (* compute angular frequency and grid dimensions *) Print["Beginning fd..."]; wo = 2 Pi fo; nx = Floor[xmax/dx + 1]; ny = Floor[h/dx + 1]; (* compute Q vector and Y admittances *) q = Table[-4 rho wo^2 dx//N,{ny}]; y = Table[(1/(s/(I wo) + beta + I wo m))//N,{x,0,xmax,dx}]; (* set up A matrices *) Print["A matrices..."]; a = Table[Table[Switch[i-j,-1,-1,0,4,1,-1,_,0], {i,ny},{j,ny}], {nx}]; Do[a[[k,1,2]]=-2; a[[k,ny,ny-1]]=-2; a[[k,1,1]]= (4 + 4 I wo rho y[[k]] dx)//N; ,{k,nx}]; (* compute B matrices *) Print["B matrices..."]; b = Table[0,{nx}]; b[[1]] = 2 Inverse[a[[1]]]; Do[ b[[k]] = Inverse[a[[k]] - b[[k-1]]],{k,2,nx-1}]; b[[nx]] = Inverse[a[[nx]] - 2 b[[nx-1]]]; (* compute C vectors *) Print["C matrices..."]; c = Table[0,{nx}]; c[[1]] = 1/2 (b[[1]] . q); Do[ c[[k]] = b[[k]] . c[[k-1]],{k,2,nx-1}]; c[[nx]] = 2 b[[nx]] . c[[nx-1]]; (* compute P vectors *) Print["P matrices..."]; p = Table[0,{nx}]; p[[nx]] = c[[nx]]//N; Do[ p[[k]] = (c[[k]] + b[[k]] . p[[k+1]])//N,{k,nx-1,1,-1}]; (* compute membrane displacement and non-differential pressure *) disp = Transpose[p][[1]] y/(I wo)//N; press = Reverse[Transpose[p/2]]; (* p is diff pressure, divide by 2 *) Print["Done fd."];]; (* :[font = input; preserveAspect; startGroup; Cclosed; nowordwrap; ] (* evaluate and plot Finite Difference method *) :[font = input; preserveAspect; startGroup; nowordwrap; ] Clear[fddisp, fdpress]; fd[rho, xmax, h, d, s, beta, m, fo, dx, fddisp, fdpress]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] A matrices... B matrices... C matrices... P matrices... :[font = input; preserveAspect; nowordwrap; ] takexmin=15; takexmax=19; p1 = surfaceplot[db[fdpress],AmbientLight->GrayLevel[0.2], LightSources->{{{1,-1,1},GrayLevel[1]}}, AxesLabel->{"x (mm)","y (mm)","|P| (dB)"}]; Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; nowordwrap; ] listplot[db[fddisp],PlotRange->{-80,20}, FrameLabel->{"Distance from Stapes (mm)","Gain (dB)"}]; listplot[phase[fddisp],FrameLabel->{"Distance from Stapes (mm)","Phase (rad)"}]; :[font = input; preserveAspect; nowordwrap; ] Show[ listplot[Re[fddisp Exp[I Pi]],AspectRatio->1/15,PlotRange->{{0,20},All}, Frame->False,PlotStyle->Thickness[.002],nodisplay], listplot[Abs[fddisp],AspectRatio->1/15,PlotRange->{{0,20},All}, Frame->False,PlotStyle->dash,nodisplay], listplot[-Abs[fddisp],AspectRatio->1/15,PlotRange->{{0,20},All}, Frame->False,PlotStyle->dash,nodisplay],display]; :[font = input; preserveAspect; nowordwrap; ] listplot[Re[fddisp Exp[I Pi]],AspectRatio->1/15,PlotRange->{{0,20},All}, Frame->False,PlotStyle->Thickness[.002]]; :[font = input; preserveAspect; nowordwrap; ] takexmax = 20; fdpresstake =Transpose[ Take[Transpose[-fdpress],pointdens*takexmax+1]]; ListDensityPlot[Re[fdpresstake],AspectRatio->1/8,Mesh->False,PlotRange->{-630000,450000}, Frame->False]; :[font = input; preserveAspect; nowordwrap; ] takexmax = 20; fdpresstake =Transpose[ Take[Transpose[-fdpress],pointdens*takexmax+1]]; ListDensityPlot[Re[fdpresstake],AspectRatio->1/10,Mesh->False,PlotRange->{-450000,450000}, FrameTicks->{Table[{pointdens*i,i},{i,0,takexmax}], Table[{i,i*h/(pointdens+1)//N},{i,0,(pointdens+1) h, (pointdens+1) h/5}]}]; :[font = input; preserveAspect; nowordwrap; ] gridplot[Re[fdpress]/Abs[fdpress],PlotRange->{-1,1}]; :[font = input; preserveAspect; endGroup; nowordwrap; ] :[font = input; preserveAspect; startGroup; nowordwrap; ] (* DEFINE: LG method *) :[font = input; initialization; preserveAspect; endGroup; nowordwrap; ] *) Clear[lg]; lg[ rho_, xmax_, h_, d_, s_, (* inputs, named as above *) beta_, m_, fo_, maxdeltak_, (* inputs, named as above *) maxdeltax_, firstguess_, (* inputs, named as above *) xlist_, (* x values, may be input or output *) k_, (* wavenumber k *) dkdx_, (* x derivative of wavenumber *) disp_, (* membrane displacement *) press_, (* fluid pressure (2-D list) *) rle_ (* relative laplace error *) ] := Module[{ (* LOCAL VARIABLES *) a, (* displacement amplitude factor *) intkdx, (* integral of k dx *) rhs, (* RHS of dispersion relation *) drhsdx, (* x-derivative of RHS of dispersion relation *) guess, (* current guess for root finding *) dum, (* dummy variable *) kval, (* current value for wavenumber k *) xval, (* current value for position x *) xstep, (* step size for x *) pressprecompute, (* precomputed value for pressure *) tanhkh, (* often-used quantity Tanh[k h] *) coshkh, (* often-used quantity Cosh[k h] *) sinhhkh, (* often-used quantity Sinh[k h] *) wo (* angular frequency in radians/s *) }, (* compute angular frequency *) Print["Beginning lg..."]; wo = 2 Pi fo; Print["solving for wavenumber..."]; (* solve for wavenumber k *) rhs = 2 rho wo^2/(s + I beta wo - m wo^2); drhsdx = D[rhs,x]; guess = firstguess; If [Length[xlist]==0, (* no previous x values, must adaptively determine them *) Print["generating x values adaptively..."]; xval=0; xstep=0; kval=0; While[xvalxmax,xval=xmax]; rhsval = rhs/.x->xval//N; If [Re[kval]>3.5, (* use Steele and Miller 1980 short-wave trick *) kval = rhsval; dkdxval = drhsdx/.x->xval, (* else evaluate in full *) kval = dum/.FindRoot[dum Tanh[dum h]==rhsval,{dum,guess}]//N; dkdxval = (drhsdx/.x->xval)/ (kval h + Tanh[kval h] - kval h Tanh[kval h]^2)//N ]; (* end if *) xstep = maxdeltak/Abs[dkdxval]//N; If[xstep>maxdeltax,xstep=maxdeltax]; guess = kval Exp [dkdxval/kval xstep]//N; If[xval == 0, xlist = {0}; k = {kval}; dkdx ={dkdxval}, AppendTo[xlist,xval]; AppendTo[k,kval]; AppendTo[dkdx,dkdxval]]; (* Print["x = ",xval," k = ",kval]; *) ], (* use previous x values *) Print["using supplied x values..."]; Do[ xval = xlist[[i]]; rhsval = rhs/.x->xval//N; kval = dum/.FindRoot[dum Tanh[dum h]==rhsval,{dum,guess}]//N; dkdxval = (D[rhs,x]/.x->xval)/ (kval h + Tanh[kval h] - kval h Tanh[kval h]^2)//N; guess = kval; If[xval == 0, k = {kval}; dkdx ={dkdxval}, AppendTo[k,kval]; AppendTo[dkdx,dkdxval]]; (* Print["x = ",xval," k = ",kval]; *) ,{i,Length[xlist]}] ]; (* end If *) Print["integrating dx..."]; (* integrate k dx, and compute some often-used quantities *) intkdx = expintlist[xlist, k]; kh = k h; tanhkh = Tanh[kh]; ko = k[[1]]; dkdxo = dkdx[[1]]; Print["computing membrane displacement..."]; (* compute membrane displacement *) a = continuous[k tanhkh/Sqrt[tanhkh + kh/(Cosh[kh]^2)]]; (* disp = I k[[1]] h a/a[[1]] *Exp[-I intkdx]; *) (* working version *) normalizer = 1/Sqrt[Tanh[ko h] + ko h/(Cosh[ko h]^2)]/(ko^2)* (dkdxo ko h + Tanh[ko h]*(-dkdxo - I ko^2 - 2 ko h dkdxo(1 - ko h Tanh[ko h])/(2 ko h + Sinh[2 ko h]) - ko h dkdxo Tanh[ko h])); disp = a h /normalizer *Exp[-I intkdx]; Print["computing pressure..."]; (* compute fluid pressure at y=0 and y=h *) pressprecompute = rho wo^2 disp/(k Sinh[kh])//N; press = Table[pressprecompute Cosh[k y]//N,{y,0,h,h}]; Print["computing relative laplace error..."]; (* compute relative laplace error at y=0 and y=h *) rle = Table[-I dkdx/(k^2) (1 - 4 kh (1 - kh tanhkh)/(2 kh + Sinh[2 kh]) + 2 k (y Tanh[k y]- h tanhkh))//N, {y,0,h,h}]; Print["done lg."]; ]; (* end Module *) (* :[font = input; preserveAspect; startGroup; nowordwrap; ] (* evaluate and plot LG method *) :[font = input; preserveAspect; startGroup; nowordwrap; ] Clear[disp1,press1,k1,dk1dx,rle1,x1]; lg[rho, xmax, h, d, s, beta, m, fo, maxdeltak, maxdeltax, guess1, x1, k1, dk1dx, disp1, press1, rle1]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] Beginning lg... solving for wavenumber... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done lg. :[font = input; preserveAspect; nowordwrap; ] xtable = x1; p1 = listplot[db[disp1],nodisplay]; xtable = Table[x,{x,0,xmax,dx}];(* table of x values *) p2 = listplot[db[fddisp],nodisplay]; Show[p1,p2, display,AspectRatio->0.7,PlotRange->{{-1,36},{-80,30}}]; :[font = input; preserveAspect; nowordwrap; ] argpointplot[k1]; :[font = input; preserveAspect; nowordwrap; ] xtable = x1; p1 = listplot[phase[disp1]/(2 Pi),nodisplay]; xtable = Table[x,{x,0,xmax,dx}];(* table of x values *) p2 = listplot[phase[fddisp]/(2 Pi),nodisplay]; Show[p1,p2, display,AspectRatio->0.7,PlotRange->{{-1,36},{-4,0.5}}]; :[font = input; preserveAspect; startGroup; nowordwrap; ] x2 = x1; (* use same x values *) Clear[disp2,press2,k2,dk2dx,rle2]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess2, x2, k2, dk2dx, disp2, press2, rle2]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] beginning adaptivewkb calculation... solving for wavenumber... using supplied x values... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done adaptivewkb. :[font = input; preserveAspect; startGroup; nowordwrap; ] nx = Length[x1]; k1h = k1 h; k2h = k2 h; tanhk1h = Tanh[k1 h]; tanhk2h = Tanh[k2 h]; Print["finding y-integrals..."]; (* find the y-integrals from closed-form approximations *) intlaplace1 = -I dk1dx/k1 press1[[2]] * (2 k1h - tanhk1h(2 + k1^2 rle1[[1]]/(I dk1dx))); intlaplace2 = -I dk2dx/k2 press2[[2]] * (2 k2h - tanhk2h(2 + k2^2 rle2[[1]]/(I dk2dx))); intp2prime = press2[[2]]/k2^2*(dk2dx k2h + tanhk2h*(-dk2dx - I k2^2 - 2 k2h dk2dx(1 - k2h tanhk2h)/(2 k2h + Sinh[2 k2h]) - k2h dk2dx tanhk2h)); intp2 = tanhk2h/k2 press2[[2]]; p = (2 intp2prime)/intp2; q = intlaplace2/intp2; r = -intlaplace1/intp2; :[font = input; preserveAspect; endGroup; nowordwrap; ] xtable=x1; :[font = input; preserveAspect; nowordwrap; ] $DefaultFont = {"Helvetica",10}; :[font = input; preserveAspect; nowordwrap; ] p1 = Show[ listplot[db[p],nodisplay,PlotStyle->dash], listplot[db[q],nodisplay,PlotStyle->gray3], listplot[db[r],nodisplay], Graphics[Text["R(x)",{8,200}]], Graphics[Text["P(x)",{8,40}]], Graphics[Text["Q(x)",{8,-22}]], display,PlotRange->{All,{-100,400}},AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)","Magnitude (dB)"}]; :[font = input; preserveAspect; nowordwrap; ] Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; startGroup; nowordwrap; ] xmax = 25; maxdeltak = .2; maxdeltax = 1; Clear[disp2,press2,k2,dk2dx,rle2,x2]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess2, x2, k2, dk2dx, disp2, press2, rle2]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] beginning adaptivewkb calculation... solving for wavenumber... generating x values adaptively... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done adaptivewkb. :[font = input; preserveAspect; nowordwrap; ] guess1 = .01-.01 I; (* starting point for root near 0 *) guess2 = .01-(Pi-.1) I//N; (* starting point for root near -I Pi *) guess3 = .01-(2 Pi-.1) I//N; (* starting point for root near -I Pi *) guess4 = .01-(3 Pi-.1) I//N; (* starting point for root near -I Pi *) guess5 = .01-(4 Pi-.1) I//N; (* starting point for root near -I Pi *) guess6 = .01-(5 Pi-.1) I//N; (* starting point for root near -I Pi *) :[font = input; preserveAspect; startGroup; nowordwrap; ] xmax = 25; maxdeltak = .2; maxdeltax = 1; Clear[disp3,press3,k3,dk3dx,rle3,x3]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess3, x3, k3, dk3dx, disp3, press3, rle3]; :[font = input; preserveAspect; startGroup; nowordwrap; ] xmax = 25; maxdeltak = .2; maxdeltax = 1; Clear[disp4,press4,k4,dk4dx,rle4,x4]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess4, x4, k4, dk4dx, disp4, press4, rle4]; :[font = input; preserveAspect; startGroup; nowordwrap; ] xmax = 25; maxdeltak = .2; maxdeltax = 1; Clear[disp5,press5,k5,dk5dx,rle5,x5]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess5, x5, k5, dk5dx, disp5, press5, rle5]; :[font = print; inactive; preserveAspect; endGroup; endGroup; endGroup; nowordwrap; ] beginning adaptivewkb calculation... solving for wavenumber... generating x values adaptively... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done adaptivewkb. :[font = input; preserveAspect; startGroup; nowordwrap; ] xmax = 25; maxdeltak = .2; maxdeltax = 1; Clear[disp6,press6,k6,dk6dx,rle6,x6]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess6, x6, k6, dk6dx, disp6, press6, rle6]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] beginning adaptivewkb calculation... solving for wavenumber... generating x values adaptively... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done adaptivewkb. :[font = input; preserveAspect; nowordwrap; ] argplot[grid_, options___] := ListPlot[Transpose[{Re[grid],Im[grid]}],options,PlotJoined->True, AspectRatio->Automatic,PlotRange->All,Frame->True,Axes->False]; argpointplot[grid_, options___] := Show[ ListPlot[Transpose[{Re[grid],Im[grid]}],PlotJoined->True,PlotRange->All, nodisplay], ListPlot[Transpose[{Re[grid],Im[grid]}],PlotRange->All,nodisplay, PlotStyle->PointSize[0.015]],options,display, AspectRatio->Automatic,PlotRange->All,Frame->True,Axes->False]; :[font = input; preserveAspect; nowordwrap; ] Show[argplot[k1,nodisplay],argplot[k2,nodisplay], argplot[k3,nodisplay],argplot[k4,nodisplay], argplot[k5,nodisplay],argplot[k6,nodisplay], display]; :[font = input; preserveAspect; startGroup; nowordwrap; ] fo = 400//N; (* frequency of input in Hz *) wo = 2 Pi fo; (* angular frequency in radians/s *) xmax = 45; maxdeltak = .2; maxdeltax = 1; Clear[disp11,press11,k11,dk11dx,rle11,x11]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess1, x11, k11, dk11dx, disp11, press11, rle11]; :[font = input; preserveAspect; startGroup; nowordwrap; ] Clear[disp11,press11,k12,dk11dx,rle11,x12]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess1, x12, k12, dk11dx, disp11, press11, rle11]; :[font = print; inactive; preserveAspect; endGroup; endGroup; nowordwrap; ] beginning adaptivewkb calculation... solving for wavenumber... generating x values adaptively... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done adaptivewkb. :[font = input; preserveAspect; nowordwrap; ] Show[argplot[k11,nodisplay],argplot[k2,nodisplay](*, argplot[k3,nodisplay],argplot[k4,nodisplay], argplot[k5,nodisplay],argplot[k6,nodisplay]*), display]; :[font = input; preserveAspect; startGroup; nowordwrap; ] x2 = x1; (* use same x values *) Clear[disp2,press2,k2,dk2dx,rle2]; adaptivewkb[rho, xmax, h, d, s, beta, m, wo, maxdeltak, maxdeltax, guess2, x2, k2, dk2dx, disp2, press2, rle2]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] beginning adaptivewkb calculation... solving for wavenumber... using supplied x values... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done adaptivewkb. :[font = input; preserveAspect; nowordwrap; ] xtable = x1; Show[ listplot[db[press1[[2]]],nodisplay], listplot[db[press1[[1]]],nodisplay], display,AspectRatio->0.8,PlotRange->{All,{0,130}}]; :[font = input; preserveAspect; nowordwrap; ] xtable = x1; Show[ listplot[db[disp2],nodisplay], listplot[db[disp2],nodisplay,PlotJoined->False,PlotStyle->PointSize[0.015]], display,AspectRatio->0.8,PlotRange->All]; Show[ listplot[phase[disp2],nodisplay], listplot[phase[disp2],nodisplay,PlotJoined->False,PlotStyle->PointSize[0.015]], display,AspectRatio->0.8,PlotRange->All]; :[font = input; preserveAspect; nowordwrap; ] xtable = x1; Show[ listplot[db[disp1],nodisplay], listplot[db[disp1],nodisplay,PlotJoined->False,PlotStyle->PointSize[0.015]], display,AspectRatio->0.8,PlotRange->{All,{-80,20}}]; Show[ listplot[phase[disp1],nodisplay], listplot[phase[disp1],nodisplay,PlotJoined->False,PlotStyle->PointSize[0.015]], display,AspectRatio->0.8,PlotRange->All]; :[font = input; preserveAspect; endGroup; nowordwrap; ] :[font = input; preserveAspect; startGroup; Cclosed; nowordwrap; ] (* DEFINE: Mode-Coupling LG method *) :[font = input; initialization; preserveAspect; endGroup; nowordwrap; ] *) mclg[ (* INPUTS TO mclg *) rho_, xmax_, h_, fo_, (* physical parameters *) dxuniform_, (* physical parameters *) x1_, (* x lists from LG solutions *) k1_, k2_, (* wavenumbers from LG solutions *) dk1dx_, dk2dx_, (* x-derivative of wavenumbers *) press1_, press2_, (* pressure solutions *) rle1_, rle2_, (* relative laplace errors *) (* OUTPUTS FROM mclg *) press_, (* fluid pressure *) disp_, (* membrane displacement *) c_ (* mixing function *) ] := Module[{ (* LOCAL VARIABLES *) intlaplace1, (* integral of L(P1) dy *) intlaplace2, (* integral of L(P2) dy *) intp2prime, (* integral of (P2)' dy *) intp2, (* integral of (P2) dy *) p, q, r, (* generic coeffs of ODE *) pgrid, qgrid, rgrid,(* generic coeffs of ODE *) xgrid, (* x values on uniform grid *) nx, nxg, (* list lengths *) wo (* angular frequency in radians/s *) }, Print["beginning mode-coupling LG calculation..."]; wo = 2 Pi fo; nx = Length[x1]; k1h = k1 h; k2h = k2 h; tanhk1h = Tanh[k1 h]; tanhk2h = Tanh[k2 h]; Print["finding y-integrals..."]; (* find the y-integrals from closed-form approximations *) intlaplace1 = -I dk1dx/k1 press1[[2]] * (2 k1h - tanhk1h(2 + k1^2 rle1[[1]]/(I dk1dx)))//N; intlaplace2 = -I dk2dx/k2 press2[[2]] * (2 k2h - tanhk2h(2 + k2^2 rle2[[1]]/(I dk2dx)))//N; intp2prime = press2[[2]]/k2^2*(dk2dx k2h + tanhk2h*(-dk2dx - I k2^2 - 2 k2h dk2dx(1 - k2h tanhk2h)/(2 k2h + Sinh[2 k2h]) - k2h dk2dx tanhk2h))//N; intp2 = tanhk2h/k2 press2[[2]]//N; p = (2 intp2prime)/intp2//N; q = intlaplace2/intp2//N; r = -intlaplace1/intp2//N; (* interpolate p,q,r onto uniform grid *) Print["interpolating onto uniform grid..."]; xgrid = Table[x,{x,0,xmax,dx}]; nxg = Length[xgrid]; pgrid = Table[expinterpolation3[x1,p,xgrid[[n]]],{n,nxg}]; qgrid = Table[expinterpolation3[x1,q,xgrid[[n]]],{n,nxg}]; rgrid = Table[expinterpolation3[x1,r,xgrid[[n]]],{n,nxg}]; (* build tridiagonal matrix to solve *) Print["building matrix..."]; lowerdiagonal= Table[1/(dx^2) - pgrid[[n+1]]/(2 dx),{n,1,nxg-1}]//N; lowerdiagonal[[nxg-1]] = 2/(dx^2)//N; upperdiagonal= Table[1/(dx^2) + pgrid[[n]]/(2 dx),{n,1,nxg-1}]//N; upperdiagonal[[1]] = 0; maindiagonal= Table[-2/(dx^2) + qgrid[[n]],{n,1,nxg}]//N; maindiagonal[[1]] = 1; rgrid[[1]]=0; (* Left-end boundary condition for ODE *) (* solve it *) Print["solving matrix..."]; cgrid = TridiagonalSolve[lowerdiagonal,maindiagonal,upperdiagonal,rgrid]; (* interpolate back to non-uniform x1 grid *) Print["interpolate back to x1 grid..."]; c = Table[expinterpolation3[xgrid,cgrid,x1[[n]]],{n,nx}]; (* construct composite solution *) press = press1 + Transpose[c Transpose[press2]]; disp = press[[2]] k1 tanhk1h/(rho wo^2)//N; Print["done."]; ]; (* END of mclg *) (* :[font = input; preserveAspect; startGroup; Cclosed; nowordwrap; ] (* evaluate and plot Mode-Coupling LG method *) :[font = input; preserveAspect; startGroup; nowordwrap; ] Clear[disp1,press1,k1,dk1dx,rle1,x1]; lg[rho, xmax, h, d, s, beta, m, fo, maxdeltak, maxdeltax, guess1, x1, k1, dk1dx, disp1, press1, rle1]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] Beginning lg... solving for wavenumber... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done lg. :[font = input; preserveAspect; nowordwrap; ] xtable = x1; listplot[db[c],PlotRange->All]; :[font = input; preserveAspect; nowordwrap; ] listplot[db[disp1],nodisplay,PlotJoined->False,PlotStyle->PointSize[0.015]], display,AspectRatio->0.8,PlotRange->{All,{-80,20}}]; Show[ listplot[phase[disp1],nodisplay], listplot[phase[disp1],nodisplay,PlotJoined->False,PlotStyle->PointSize[0.015]], display,AspectRatio->0.8,PlotRange->All]; :[font = input; preserveAspect; startGroup; nowordwrap; ] Clear[disp2,press2,k2,x2,dk2dx,rle2]; lg[rho, xmax, h, d, s, beta, m, fo, maxdeltak, maxdeltax, guess2, x2, k2, dk2dx,disp2, press2, rle2]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] Beginning lg... solving for wavenumber... integrating dx... computing membrane displacement... computing pressure... computing relative laplace error... done lg. :[font = input; preserveAspect; startGroup; nowordwrap; ] (* run adaptive mode-coupling routine with direct solve for c(x) *) dxuniform = 1/11; Clear[disp3,press3,c]; mclg[rho, xmax, h, fo, dxuniform, x1, x2, k1, k2, dk1dx, dk2dx, press1, press2,rle1, rle2, press3, disp3, c]; :[font = message; inactive; preserveAspect; endGroup; nowordwrap; ] 1 Power::infy: Infinite expression - encountered. 0 :[font = input; preserveAspect; startGroup; Cclosed; nowordwrap; ] (* saving and reading partial results to a file *) :[font = input; preserveAspect; nowordwrap; ] (* save all computation results to archive file *) Timing[Save["math/steeleneely/archivestate5", rho, xmax, h, d, s, beta, m, wo, pointdens, dx, nx, ny, guess2, guess1, xtable, fddisp, fdpress, x1, k1, dk1dx, disp1, press1, rle1, x2, k2, dk2dx, disp2, press2, rle2];][[1]] :[font = input; preserveAspect; endGroup; nowordwrap; ] (* load all computation results from archive file *) Timing[<<"math/steeleneely/archivestate5";][[1]] :[font = input; preserveAspect; endGroup; nowordwrap; ] (* compare finite diff, wkb and MCwkb methods *) Show[listplot[db[disp1],PlotStyle->dash,nodisplay], listplot[db[disp3],PlotStyle->gray3,nodisplay], (* listplot[db[fddisp],nodisplay], *) PlotRange->{All,{-80,20}},display, FrameLabel->{"Distance from stapes (mm)","Gain (dB)"}]; Show[listplot[phase[disp1],PlotStyle->dash,nodisplay], listplot[phase[disp3],PlotStyle->gray3,nodisplay], (* listplot[phase[fddisp],nodisplay], *) PlotRange->All,display, FrameLabel->{"Distance from stapes (mm)","Phase (rad)"}]; :[font = input; preserveAspect; startGroup; nowordwrap; ] (* multi-frequency comparisons *) :[font = input; preserveAspect; startGroup; nowordwrap; ] (* PHYSICAL PARAMETERS *) rho = .001; (* density of water in g/mm3 *) h = 1; (* height of scalae in mm *) d = 5; (* characteristic length of cochlea in mm *) s = 10 10^6 Exp [-x/d]; (* membrane stiffness/area in g/(s2mm2) *) beta = 2; (* membrane damping/area in g/smm2 *) m = 1.5 10^(-3); (* membrane mass/area in g/mm2 *) (* COMPUTATIONAL PARAMETERS *) maxdeltak = .1; (* maximum allowed change in k per step *) maxdeltax = .5; (* maximum step size *) guess1 = .01-.01 I; (* starting point for root near 0 *) guess2 = .01-(Pi-.1) I//N; (* starting point for root near -I Pi *) (* FREQUENCY SWEEP *) fstart = 400; (* starting frequency in Hz *) findexinit = 8; (* initial findex *) numf = 10; (* number of frequencies at Sqrt[2] steps *) xmaxarray = {35,35,35,32,27,23,20,15,12,8}; (* length of cochlea to consider for each freq *) dxuniformarray = {1/11,1/11,1/11,1/11,1/11,1/11,1/11,1/20,1/20,1/20}; (* uniform grid spacing for mcwkb soln *) fddxarray = {1/9,1/9,1/9,1/9,1/9,1/11,1/11,1/11,1/13,1/13}; (* point spacing for finite diff soln *) fdxarray = Table[0,{numf}]; fddisparray = Table[0,{numf}]; disp1array = Table[0,{numf}]; disp2array = Table[0,{numf}]; disp3array = Table[0,{numf}]; x1array = Table[0,{numf}]; carray = Table[0,{numf}]; Do[ dx = fddxarray[[findex]]; xmax = xmaxarray[[findex]]; nx = Floor[xmax/dx + 1]; (* number of grid points in x direction (M) *) ny = Floor[h/dx + 1]; fo = fstart Sqrt[2]^(findex-1)//N; pointdens = Floor[1/dx]; (* point density *) dxuniform = dxuniformarray[[findex]]; Print["***********************"]; Print["Beginning new frequency\n"]; Print["fo = ",fo]; Print["xmax = ",xmax]; Print["dx = ",dx]; Print["dxuniform = ",dxuniform]; Print[" ",dxuniform]; Run["date"]; Print[""]; (* compute finite difference solution *) Clear[fddisp, fdpress]; fd[rho, xmax, h, d, s, beta, m, fo, dx, fddisp, fdpress]; fdx = Table[x,{x,0,xmax,dx}]; Run["date"]; Print[""]; Clear[disp1,press1,k1,dk1dx,rle1,x1]; lg[rho, xmax, h, d, s, beta, m, fo, maxdeltak, maxdeltax, guess1, x1, k1, dk1dx,disp1, press1, rle1]; Run["date"]; Print[""]; x2 = x1; Clear[disp2,press2,k2,dk2dx,rle2]; lg[rho, xmax, h, d, s, beta, m, fo, maxdeltak, maxdeltax, guess2, x2, k2, dk2dx, disp2, press2, rle2]; Run["date"]; Print[""]; (* run adaptive mode-coupling routine *) Clear[disp3,press3,c]; mclg[rho, xmax, h, fo, dxuniform, x1, k1, k2, dk1dx, dk2dx, press1, press2,rle1, rle2, press3, disp3, c]; Run["date"]; Print[""]; (* store results *) fdxarray[[findex]] = fdx; fddisparray[[findex]] = fddisp; disp1array[[findex]] = disp1; disp2array[[findex]] = disp2; disp3array[[findex]] = disp3; x1array[[findex]] = x1; carray[[findex]] = c; (* plot magnitudes *) Show[ ListPlot[Transpose[{x1array[[findex]],db[disp3array[[findex]]]}], PlotJoined->True,PlotStyle->gray3,nodisplay], ListPlot[Transpose[{x1array[[findex]],db[disp1array[[findex]]]}], PlotJoined->True,PlotStyle->dash,nodisplay], ListPlot[Transpose[{fdxarray[[findex]],db[fddisparray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay], PlotRange->{All,{-80,20}},display,Frame->True,Axes->False, AspectRatio->0.8, FrameLabel->{"Distance from stapes (mm)","Gain (dB)"}]; ,{findex,findexinit,numf}]; :[font = input; preserveAspect; startGroup; nowordwrap; ] Print["Saving data to a file..."]; Save["/ufs/physcmp/lloyd/math/mcwkb/multif", numf, fdxarray, fddisparray, disp1array, disp2array, disp3array, x1array, carray]; Print["done."]; :[font = print; inactive; preserveAspect; endGroup; nowordwrap; ] Saving data to a file... done. :[font = input; preserveAspect; nowordwrap; ] <True,PlotStyle->dash,nodisplay], ListPlot[Transpose[{x1array[[8]],phase[disp2array[[8]]]/(2 Pi)}], PlotJoined->True,PlotStyle->gray4,nodisplay], ListPlot[Transpose[{x1array[[8]],phase[disp3array[[8]]]/(2 Pi)}], PlotJoined->True,PlotStyle->solid,nodisplay], Graphics[Text["Composite",{6,-2},{1,0}]], Graphics[Text["Traveling-Wave",{6,-2.5},{1,0}]], Graphics[Text["Cut-Off",{6,-3},{1,0}]], Graphics[{solid,Line[{{6.5,-2},{8,-2}}]}], Graphics[{dash,Line[{{6.5,-2.5},{8,-2.5}}]}], Graphics[{gray4,Line[{{6.5,-3},{8,-3}}]}], PlotRange->{{-.5,15.5},All},display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Phase of\nDisplacement Ratio (cycles)"}]; Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; startGroup; nowordwrap; ] 3.2*Sqrt[2]//N :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 4.525483399593905 ;[o] 4.52548 :[font = input; preserveAspect; nowordwrap; ] p1 = Show[ ListPlot[Transpose[{x1array[[8]],db[disp3array[[8]]]}], PlotJoined->True,PlotStyle->solid,nodisplay], PlotRange->{{-.5,15.5},{-81,30}},display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Magnitude of\nDisplacement Ratio (dB)"}]; Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; nowordwrap; ] p1 = Show[ ListPlot[Transpose[{x1array[[6]],db[carray[[6]]*disp2array[[6]]]}], PlotJoined->True,PlotStyle->gray4,nodisplay], ListPlot[Transpose[{x1array[[6]],db[disp1array[[6]]]}], PlotJoined->True,PlotStyle->dash,nodisplay], ListPlot[Transpose[{fdxarray[[6]],db[fddisparray[[6]]]}], PlotJoined->True,PlotStyle->solid,nodisplay], Graphics[Text["composite",{22,12},{-1,0}]], Graphics[Text["traveling wave",{22,0},{-1,0}]], Graphics[Text["cut-off",{22,-12},{-1,0}]], Graphics[{solid,Line[{{19,12},{21.3,12}}]}], Graphics[{dash,Line[{{19,0},{21.3,0}}]}], Graphics[{gray4,Line[{{19,-12},{21.3,-12}}]}], PlotRange->{{-.5,30.5},{-81,30}},display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Magnitude of\nDisplacement Ratio (dB)"}]; Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; nowordwrap; ] numf=10; :[font = input; preserveAspect; nowordwrap; ] fddisparray[[9]] = fddisp; fdxarray[[9]] = xtable; :[font = input; preserveAspect; nowordwrap; ] fddisparray[[10]] = fddisp; fdxarray[[10]] = xtable; :[font = input; preserveAspect; nowordwrap; ] disp3array[[10]] = disp3; x1array[[10]] = x1; carray[[10]] = c; :[font = input; preserveAspect; startGroup; nowordwrap; ] p1 = Show[ Table[ListPlot[Transpose[{x1array[[findex]],db[disp1array[[findex]]]}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],db[fddisparray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], Graphics[Text["0.4",{31,16}]], Graphics[Text["0.8",{25,17}]], Graphics[Text["1.6",{19,18}]], Graphics[Text["3.2",{12.5,19}]], Graphics[Text["6.4 kHz",{6.5,20}]], PlotRange->{All,{-80,30}},display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Magnitude of\nDisplacement Ratio (dB)"}]; :[font = input; preserveAspect; nowordwrap; ] p2 = Show[ Table[ListPlot[Transpose[{x1array[[findex]],phase[disp1array[[findex]]]/(2 Pi)}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],phase[fddisparray[[findex]]]/(2 Pi)}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], Graphics[Text["0.4",{33.4,-1.3}]], Graphics[Text["0.8",{27.5,-1.7}]], Graphics[Text["1.6",{20.6,-1.8}]], Graphics[Text["3.2",{14.1,-1.85}]], Graphics[Text["6.4",{7.8,-1.9}]], Graphics[Text["kHz",{8,-2.1}]], PlotRange->All,display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Phase of\nDisplacement Ratio (cycles)"}]; :[font = input; preserveAspect; startGroup; nowordwrap; ] Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; endGroup; nowordwrap; ] Show[Graphics[Rectangle[{0,0},{1,1},p2]]]; :[font = input; preserveAspect; endGroup; nowordwrap; ] single-mode wkb :[font = input; preserveAspect; startGroup; nowordwrap; ] p3 = Show[ Table[ListPlot[Transpose[{x1array[[findex]],db[disp3array[[findex]]]}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],db[fddisparray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], Graphics[Text["0.4",{31,16}]], Graphics[Text["0.8",{25,17}]], Graphics[Text["1.6",{19,18}]], Graphics[Text["3.2",{12.5,19}]], Graphics[Text["6.4 kHz",{6.5,20}]], PlotRange->{All,{-80,30}},display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Magnitude of\nDisplacement Ratio (dB)"}]; :[font = input; preserveAspect; nowordwrap; ] p4 = Show[ Table[ListPlot[Transpose[{x1array[[findex]],phase[disp3array[[findex]]]/(2 Pi)}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],phase[fddisparray[[findex]]]/(2 Pi)}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], Graphics[Text["0.4",{33.4,-1.3}]], Graphics[Text["0.8",{27.5,-1.7}]], Graphics[Text["1.6",{20.6,-1.8}]], Graphics[Text["3.2",{14.1,-1.85}]], Graphics[Text["6.4",{7.8,-1.9}]], Graphics[Text["kHz",{8,-2.1}]], PlotRange->All,display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Phase of\nDisplacement Ratio (cycles)"}]; :[font = input; preserveAspect; endGroup; nowordwrap; ] Show[Graphics[Rectangle[{0,0},{1,1},p3]]]; Show[Graphics[Rectangle[{0,0},{1,1},p4]]]; :[font = input; preserveAspect; startGroup; nowordwrap; ] Show[ Table[ListPlot[Transpose[{x1array[[findex]],db[disp3array[[findex]]]}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],db[fddisparray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], PlotRange->{All,{-80,20}},display,Frame->True,Axes->False, AspectRatio->0.8, FrameLabel->{"Distance from stapes (mm)","Gain (dB)"}]; Show[ Table[ListPlot[Transpose[{x1array[[findex]],phase[disp3array[[findex]]]}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],phase[fddisparray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], PlotRange->All,display,Frame->True,Axes->False, AspectRatio->0.8, FrameLabel->{"Distance from stapes (mm)","Phase (rad)"}]; :[font = input; preserveAspect; endGroup; endGroup; nowordwrap; ] mode-coupling wkb, with direct solution for c(x) :[font = input; preserveAspect; nowordwrap; ] p1 = Show[ Table[ListPlot[Transpose[{x1array[[findex]],db[disp1array[[findex]]]}], PlotJoined->True,PlotStyle->dash,nodisplay],{findex,numf}], Table[ListPlot[Transpose[{fdxarray[[findex]],db[fddisparray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], Graphics[Text["0.4",{31,16}]], Graphics[Text["0.8",{25,17}]], Graphics[Text["1.6",{19,18}]], Graphics[Text["3.2",{12.5,19}]], Graphics[Text["6.4 kHz",{6.5,20}]], PlotRange->{{-1,36},{-80,30}},display,Frame->True,Axes->False, AspectRatio->0.7, FrameLabel->{"Distance from stapes (mm)"," Magnitude of\nDisplacement Ratio (dB)"}]; :[font = input; preserveAspect; startGroup; nowordwrap; ] Length[x1array] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 10 ;[o] 10 :[font = input; preserveAspect; startGroup; nowordwrap; ] numf = 10; p1 =Show[ Table[ListPlot[Transpose[{x1array[[findex]],db[carray[[findex]]]}], PlotJoined->True,PlotStyle->solid,nodisplay],{findex,numf}], Graphics[Text["0.4",{33,820}]], Graphics[Text["0.8",{30,650}]], Graphics[Text["1.6",{24,465}]], Graphics[Text["3.2",{18,285}]], Graphics[Text["6.4 kHz",{11,105},{-1,0}]], PlotRange->{All,{-160,900}},display,Frame->True,Axes->False, FrameLabel->{"Distance from stapes (mm)", "Magnitude of c(x) (dB)"}, AspectRatio->0.8]; Show[Graphics[Rectangle[{0,0},{1,1},p1]]]; :[font = input; preserveAspect; endGroup; nowordwrap; ] mixing factor c(x) :[font = input; preserveAspect; startGroup; nowordwrap; ] (* saving and reading partial results to a file *) :[font = input; preserveAspect; nowordwrap; ] (* save all computation results to archive file *) Timing[Save["math/steeleneely/archivestate5", rho, xmax, h, d, s, beta, m, wo, pointdens, dx, nx, ny, guess2, guess1, xtable, fddisp, fdpress, k1, dk1dx, disp1, press1, rle1, k2, dk2dx, disp2, press2, rle2];][[1]] :[font = input; preserveAspect; endGroup; endGroup; nowordwrap; ] (* load all computation results from archive file *) Timing[<<"math/steeleneely/archivestate5";][[1]] :[font = input; preserveAspect; nowordwrap; ] Please logout without saving if you need to use anvil. thanks, Lloyd ^*