Hi

I've made a small script to play around with Bessel functions of the first kind... but this is very basic and I'm wondering if there's a smarter way.

Below is the script I made (the function + a small test which plots the result).

Best regards,
Claus

// bessel_test.sce

function  z=Jn(n,x)
    // The following power series approximates the nth Bessel
    // function of the first kind for each input x
    // In acoustics x = 2ka, defines the iput frequency k = omega / c and size
    // of the piston radiator
    powerseries  =  0;
    for  m=0:19  // actually it should be infinity, but 10 approximates OK ...
        powerseries_m  =  ((-1)^m  /  (factorial(m)  *  factorial(m  +  n)))  * 
 (x/2)^(2*m);
        powerseries  =  powerseries  +  powerseries_m;  // sum the powerseries
    end
    z  =  ((x/2)^n)  .*  powerseries;
endfunction

x_array  =  0:0.1:9.9;  // define 100 points on the x-axis
                     // with 2ka (x) from 0 to 10, and with m = 0-19,
                     // this approximation is reasonably good for 2ka<  10

z_0_array  =  Jn(0,x_array);  // Calculate J0

z_1_array  =  Jn(1,x_array);  // Calculate J1

scf();
a  =  gca();
plot(x_array,z_0_array,'-b');
plot(x_array,z_1_array,'-r');
xtitle("Bessel functions","x (2ka)","output (z)");
legend("J0","J1");


_______________________________________________
users mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/users

Reply via email to