/*
* cc bessel.c -lm 
* Include the first 5 terms
*/
        #include <stdio.h>
        # include <math.h>
main()   {
        double rr, drr;
        double j0(double x);
        double j1(double x);
        double j0x,j1x;
        double G, time, uu, lam;
        int i;

        G = 1./3.;

/*  Time 1  */
        time = 0.05;
        
        rr = 0.0;
        drr = 0.01;
        i=1;
        while (i <= 101) {
        lam = 2.405;
        uu = j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 5.520;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 8.654;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 11.792;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 14.931;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);

        uu = -2.*G*uu + G/4. *(1.-rr*rr);
        printf("%f %f %f\n",rr,time,uu);
        rr = rr+drr;
        i++;
        }

/*  Time 1  */
        time = 0.1;

        rr = 0.0;
        drr = 0.01;
        i=1;
        while (i <= 101) {
        lam = 2.405;
        uu = j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 5.520;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 8.654;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 11.792;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 14.931;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);

        uu = -2.*G*uu + G/4. *(1.-rr*rr);
        printf("%f %f %f\n",rr,time,uu);
        rr = rr+drr;
        i++;
        }

/*  Time 1  */
        time = 0.2;

        rr = 0.0;
        drr = 0.01;
        i=1;
        while (i <= 101) {
        lam = 2.405;
        uu = j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 5.520;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 8.654;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 11.792;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 14.931;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);

        uu = -2.*G*uu + G/4. *(1.-rr*rr);
        printf("%f %f %f\n",rr,time,uu);
        rr = rr+drr;
        i++;
        }

/*  Time 1  */
        time = 0.4;

        rr = 0.0;
        drr = 0.01;
        i=1;
        while (i <= 101) {
        lam = 2.405;
        uu = j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 5.520;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 8.654;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 11.792;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 14.931;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);

        uu = -2.*G*uu + G/4. *(1.-rr*rr);
        printf("%f %f %f\n",rr,time,uu);
        rr = rr+drr;
        i++;
        }

/*  Time 1  */
        time = 0.75;

        rr = 0.0;
        drr = 0.01;
        i=1;
        while (i <= 101) {
        lam = 2.405;
        uu = j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 5.520;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 8.654;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 11.792;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);
        lam = 14.931;
        uu = uu + j0(lam*rr)*exp(-lam*lam*time)/lam/lam/lam/j1(lam);

        uu = -2.*G*uu + G/4. *(1.-rr*rr);
        printf("%f %f %f\n",rr,time,uu);
        rr = rr+drr;
        i++;
        }

        
        }