/* Imagine a cylinder filled with a liquid, surrounded by solid 1, in turn surrounded by solid 2
   each of a set thickness.

   For example, the liquid could be drilling mud in a bore-hole, solid 1 the drilling mud invaded zone
   and solid 2 the untouched rock through which the borehole passes.
 */


/* Parameters */

a:       0.2$   /* Borehole diameter, metres */

alpha_1: 1500$  /* Sound velocity in liquid, metres/sec*/
rho_1:   1000$     /* Fluid density, kg/m^3*/


alpha_2: 1830$  /* Compressional velocity in solid 1, meters/sec */
beta_2:  914$   /* Shear velocity in solid 1, metres/sec */
rho_2:   2000$     /* Density of solid 1, kg/m^3 */

/* freq:    2500$ */  /* Excitation frequency, Hz */
/* k:        300$ */  /* Wavenumber */

omega:   2 * %pi * freq$

v1:      (k^2 - omega^2/(alpha_1^2))^(1/2)$=20
v2:      (k^2 - omega^2/(alpha_2^2))^(1/2)$
v3:      (k^2 - omega^2/(beta_2^2))^(1/2)$

h:       float ( rho_1/rho_2 * ( a * omega / beta_2 )^2 )$

det111: -2 * a * v2 * bessel_k( n+1, a*v2 ) - bessel_k ( n, a * v2 ) * (2 * n * (n-1) + a * a * (2*k*k - omega*omega / (beta_2*beta_2) ))$
det112:  2 * n * (-a * v3 * bessel_k ( n+1, a * v3 ) + (n-1) * bessel_k ( n, a * v3))$
det113:  2 * k * a * (a *v3 * bessel_k ( n+1, a * v3 ) + bessel_k ( n, a * v3 ) * (n * (n-1) + (a * v3)^2))$
det121:  2 * n * (-a * v2 * bessel_k ( n+1, a * v2 ) + (n-1) * bessel_k ( n, a * v2))$
det122: -2 * a * v3 * bessel_k ( n+1, a * v3 ) - bessel_k ( n, a * v3 ) * (2 * n * (n-1) + (a * v3)^2)$
det123:  2 * a * k * n * (a * v3 * bessel_k ( n+1, a * v3 ) - (n-1) * bessel_k ( n, a * v3 ))$
det131:  2 * k * a * (-a * v2 * bessel_k ( n+1, a * v2 ) + n * bessel_k ( n, a * v2 ))$
det132: -a * k * n * bessel_k ( n, a * v3 )$
det133:  a * a * (omega*omega  / (beta_2*beta_2) - 2 * k * k) * ( - a * v3 * bessel_k ( n+1, a * v3 ) + n * bessel_k ( n, a * v3 ) )$

det211:  a * v2 * bessel_k ( n+1, a * v2 ) - n * bessel_k ( n, a * v2 )$
det212:  n * bessel_k ( n, a * v3 )$
det213:  k * a * (-a * v3 * bessel_k ( n+1, a * v3 ) + n * bessel_k ( n, a * v3 ))$
det221:  2 * n * (-a * v2 * bessel_k ( n+1, a * v2 ) + (n-1) * bessel_k ( n, a * v2 ))$
det222: -2 * a * v3 * bessel_k ( n+1, a * v3 ) - bessel_k ( n, a * v3 ) * (2 * n * (n-1) + (a * v3)^2)$
det223:  2 * a * k * n * (a * v3 * bessel_k ( n+1, a * v3 ) - (n-1) * bessel_k ( n, a * v3 ))$
det231:  2 * k * a * (-a * v2 * bessel_k ( n+1, a * v2 ) + n * bessel_k ( n, a * v2 ))$
det232: -a * k * n * bessel_k ( n, a * v3 )$
det233:  a * a * (omega*omega  / (beta_2*beta_2) - 2 * k * k) * ( a * v3 * bessel_k ( n+1, a * v3 ) + n * bessel_k ( n, a * v3 ) )$

mat1:    matrix ( [det111, det112, det113],
		  [det121, det122, det123],
		  [det131, det132, det133] )$

det1:    determinant ( mat1 )$

mat2:    matrix ( [det211, det212, det213],
		  [det221, det222, det223],
		  [det231, det232, det233] )$

det2:    determinant ( mat2 )$

deltan:   det1 * ( v1 * a * bessel_i ( n+1, v1 * a ) + n * bessel_i ( n, v1 * a ) ) + h * det2 * bessel_i ( n, v1 * a )$
deltan1: deltan,n=1;

/*block([w:0], for u:100 thru 150 do (for v:1000 thru 1004 do ( w:ev(deltan1,freq=u,k=v),print(u, " ", v, " ", float(w)))));*/

/*plot3d(deltan1,[freq,1000,1005],[k,300,350]);*/
/*plot3d(deltan1,[freq,100,2000],[k,15,50]);*/
/*plot3d(deltan1,[freq,50,70],[k,1,10]);*/
 /* plot3d(deltan1,[freq,1,112],[k,0.8,15]);*/
/*plot3d(deltan1,[freq,1,112],[k,0.8,100]);*/
/*plot3d(deltan1,[freq,1,112],[k,0.78,15]);*/

/* s1:      solve ( deltan, k ), n=1$ */