/* 
The following program is a c version of the effmass.f program .
It calculates the internal mode frequencies of a right circular
cylinder using the method outlined by J.R. Hutchinson (JR Hutchinson,
J. Appl. Mech., Vol. 47, pp 901-907, Dec. 1980.). It also calculates
the "effective mass coefficients" for each mode using  the method
developed by Aaron Gillespie.


Modification History : 

   Original Fortran code by J.R. Hutchinson

   Modified and extended to calculate effective mass coefficients
   by Aaron Gillespie. (effmass.3.2.f)

   Rewritten by Kent Blackburn in C, with modifications that eliminated
   various numerical instabilites and poor convergence problems 
   displayed by effmass.3.2.f.

   Received by John Carri from Kent Blackburn November 21, 1996.
   
   function opening_message() added Nov 22.
   Modified to read input from file or keyboard , chosen by
   command-line argument, Nov 22 1996. Saved as testmass3.c .	

   Received by Kent Blackburn from John Carri January 1997.
   changed usage of reserve inline variable to in_line (ANSI violation)
   fixed fflush(infile)/fclose(infile) bug that caused core dump w/ gcc
   added USAGE descriptor, Feb. 13 1997

   Modified July 22, 1997 by J.K.Blackburn to add length, radius, beamspot
   poisson number and young's modulus to the header of the output file
   allowing the identification of the parameters used in generating the
   results. This is to make the output compatable with the new C++ noise
   class library developed by SURF Project.
*/

#include <stdio.h>
#include <math.h>

#define PI	3.14159265358979323846264338327950

#define MIN( a, b )  ((a) < (b) ? (a) : (b))
#define MAX( a, b )  ((a) > (b) ? (a) : (b))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

static int    nz,nr,n,iterm;
static double h,sd,v,dnor;
static int    ndr,ndz;
static double shear,rad,tmass,dens,speed;
static double d[41][41],rj[51][6];

static double a[41],b[41][41],cp[41],as[41],bs[41][41],
       aba[41][41],bba[41],ah[41],bh[41][41],ch[41],aa[41],
       aas[41],bas[41],sba[41],ba[41],das[41],sda[41],
       da[41],bk[41],bp[41],bj[41],dk[41],dp[41],dj[41],
       ab[41],ab2[41],bbs[41],sbb[41],bb[41],
       dbs[41],sdb[41],db[41],ss[41],cs[41],csd[41],
       ao2[41],cba[41][41],abas[41][41],bbas[41][41];

static double veca[41][2001],vecb[41][2001],vecc[41][2001],
              norm[2001];

static double psiba[41],qjbba[41],qjpba[41],psida[41],qjpda[41],
       qjbda[41],psiab[41],qjpab[41],qjbab[41],cosdbh[41],
       cosbbh[41],cosaah[41],sndbh[41],snbbh[41],snaah[41];

static double psibar[41],qjbbar[41],qjpbar[41],psidar[41],qjpdar[41],
       qjbdar[41],psiabr[41],qjpabr[41],qjbabr[41],cosdbz[41],
       cosbbz[41],cosaaz[41],sndbz[41],snbbz[41],snaaz[41];

static int debug = 1;
static int dbgall = 0;

int main (int argc, char* argv[])
{
 int i;
 int nodes, numr, numz, numdr, numdz, numrf;
 double width, radius, spotsize, beginf, endf, deltaf, youngs,
        poisson, density, resonf[4001], alpha[4001];
 int parity[2001],order[2001];
 int read_keyboard = 1;
 FILE *infile, *outfile;
 char in_line[256];
 char infilename[256];
 const char outfilename[] = "testmass.out";
 void buildeffmassfreq();
 void opening_message(char name[]);
 
/* LIGO Testmass Default Variable Initialization */

 numrf = 0;
 numdr = 40;
 numdz = 40;
 
 nodes = 0; width = 10.0; radius = 12.5; spotsize = 2.2;
 beginf = 3000; endf = 50000; deltaf = 25.0;
 youngs = 0.717e12; poisson = 0.16; density = 2.2;

 /* request parameters from the user */ 

 opening_message(argv[0]);

 for (i=1;i<argc;i++) {
        if( strcmp(argv[i], "-file") == 0 && i < argc-1){
		strcpy(infilename,argv[++i]);
		if( (infile = fopen(infilename,"r")) == NULL){
			printf("Could not open data file \"%s\"\n",infilename);
			printf("continuing with interactive data input...\n");
		}
		read_keyboard = 0;
	}else
		read_keyboard = 1;
 }
 

 if(read_keyboard){   /* read input from keyboard, prompting for inputs...*/

 	printf(" YOU MUST SUPPLY INPUTS FOR ALL REQUESTED PARAMETERS!\n");
 	printf(" ALL PARAMETERS MUST BE GIVEN IN CGS UNITS!\n");
  
 	printf(" \nEnter the number of nodal lines on the faces of the cylinder:\n");
 	printf(" [effective mass can only be calculated for order 0] \n");
 	printf(" [allowed orders are 0,1,2,3,4]\n");
 	scanf("%d",&nodes);
 
 	printf(" \nEnter the length of the cylinder in centimeters:\n");
 	scanf("%lf",&width);
 
 	printf(" \nEnter the radius of the cylinder face in centimeters:\n");
 	scanf("%lf",&radius);
 
 	printf(" \nEnter the beam spotsize radius in centimeters:\n");
 	printf(" [radius where the field is 1/e of its maximum]\n");
 	scanf("%lf",&spotsize);
 
 	printf(" \nEnter the starting frequency for searching out resonant modes:\n");
 	printf(" [should be less than or of order 1000 depending on cylinder]\n");
 	scanf("%lf",&beginf);
 
 	printf(" \nEnter the stopping frequency for searching out resonant modes:\n");
 	printf(" [should be greater than 200K Hz to saturate thermal noise] \n");
 	printf(" [should be greater than 30K Hz to search for lowest modes] \n");
 	scanf("%lf",&endf);
 
 	printf(" \nEnter the frequency stepsize to search for resonant modes:\n");
 	printf(" [between 10 and 100 Hz works best]\n");
 	printf(" [larger is faster, but may miss some modes]\n");
 	scanf("%lf",&deltaf);
 
 	printf(" \nEnter the Youngs Modulus of the material:\n");
 	printf(" [for fused silica it is 0.717e12]\n");
 	scanf("%lf",&youngs);
 
 	printf(" \nEnter the Poisson Ratio of the material:\n");
 	printf(" [for fused silica it is 0.16]\n");
 	scanf("%lf",&poisson);
 	
 	printf(" \nEnter the density of the material:\n");
 	printf(" [for fused silica it is 2.2]\n");
 	scanf("%lf",&density);
 
 	printf(" \nEnter the number of radial terms used in the expansion:\n");
 	printf(" [must be an integer between 4 and 40]\n");
 	printf(" [smaller is faster, but WILL miss some modes]\n");
 	printf(" [20 seems to be an optimal value, I like 30]\n");
 	scanf("%d",&numr);
 
 	printf(" \nEnter the number of axial terms used in the expansion:\n");
 	printf(" [must be an integer between 4 and 40]\n");
 	printf(" [smaller is faster, but WILL miss some modes]\n");
 	printf(" [20 seems to be an optimal value, I like 30]\n");
 	scanf("%d",&numz);

  } else {  /* read from specified data file  */

	fgets(in_line,120,infile);
	sscanf(in_line,"%d",&nodes);
 
	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&width);
 
	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&radius);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&spotsize);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&beginf);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&endf);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&deltaf);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&youngs);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&poisson);

	fgets(in_line,120,infile);
	sscanf(in_line,"%lf",&density);

	fgets(in_line,120,infile);
	sscanf(in_line,"%d",&numr);

	fgets(in_line,120,infile);
	sscanf(in_line,"%d",&numz);

        fflush(infile);
        fclose(infile);
 }

 if((outfile = fopen(outfilename,"w")) == NULL){
	fprintf(stderr,"Could not open output file \"%s\". Aborting.\n",outfile);
	exit(1);
 }
 fprintf(stderr,"\n\nCalculating resonant frequencies and effective masses. ");
 fprintf(stderr,"Please wait.....\n\n");
 buildEffmassFreq(nodes,width, radius, spotsize, beginf,
        endf, deltaf, youngs, poisson, density, numr, numz,
        numdr, numdz, resonf, alpha, parity, order, &numrf);

 /* The following parameters identify results file */
 fprintf(outfile,"%f \n",width);
 fprintf(outfile,"%f \n",radius);
 fprintf(outfile,"%f \n",density);
 fprintf(outfile,"%f \n",spotsize);
 fprintf(outfile,"%f \n",poisson);
 fprintf(outfile,"%f \n",youngs);
 
 for (i=0; i<numrf; i++)
    {
     if (n == 0)
       fprintf(outfile,"%f %f  %d %d %d \n",
               resonf[i],alpha[i],nodes,parity[i],order[i]);
     else
       fprintf(outfile,"%f  %d %d %d \n",resonf[i],nodes,parity[i],order[i]);

    }
 fprintf(stderr,"\nCalculation complete. Output in file \"%s\"\n\n",outfilename);

 fflush(outfile);
 
 fclose(outfile);
 return(0);
}
 

buildEffmassFreq(NODES,L,R,S,FB,FE,D,Y,V,RHO,NR,NZ,NDR,NDZ,
                 reson,alpha,parity,order,NRF)
int NODES,NR,NZ,NDR,NDZ,*NRF;
double L,R,S,FB,FE,D,Y,V,RHO;
double *reson,*alpha;
int *parity,*order;
{
 FILE *file;
 static int nd = 40;
 static int i,j,nf,nfreq,iacc,iwrt,ifreq,ierr,kff,nff;
 static double oma,omm,del,spot,omega,om1,om2,dom1,dom2,omroot,acc;
 static double ada,adb,omc,bc,s2v,bmc,amt,y,conv;
 static double bcsav[2001],freq[2001],weight[2001];
 extern int n,iterm;
 extern double h,rad,dens,dnor,shear,speed,v,d[41][41],rj[51][6],norm[2001];
 void effmass();
 double findroot(),func();

 *NRF = 0;
 h = L / (2 * R);
 rad = R;
 spot = S;
 oma = FB; omm = FE;
 del = D;
 y = Y;
 v = V;
 dens = RHO;
 nr = NR; nz = NZ;
 ndr = NDR; ndz = NDZ;
 shear = Y / 2 / (1 + V);
 speed = sqrt(shear/dens);
 conv = 2 * PI * rad  / speed; 
 oma *= conv; omm *= conv; del *= conv;
 n = NODES;
 acc = 1.0e-24;
 dnor = 1.0;
 
 if((file = fopen("zeros.dbl","r")) == NULL){
 	fprintf(stderr,"Could not find data file \"zeros.dbl\". Aborting.\n");
	exit(1);
 }
 
 for (i=1; i<=50; i++)
    fscanf(file,"%lf %lf %lf %lf %lf\n",
          &rj[i][1],&rj[i][2],&rj[i][3],&rj[i][4],&rj[i][5]); 
 fclose(file);

 for (iterm = 0; iterm <=1; iterm++)
    {
     dnor = 1.0;
     nfreq = 0;
     om1 = omega = oma;
     dom1 = func(om1);
     dnor *= pow( 1.0/fabs(func(om1+del/2.0)),1.0/nr);
     printf("DNOR: %f\n",dnor);
     while (omega <= omm)
          {
           om2 = omega + del;
           dom2 = func(om2);
           if (dom1 * dom2 < 0.0) 
             {
              omroot = findroot(func,om1,om2,acc);
              if (fabs(dom1 = func(omroot)) < 1.0e-6) 
                {
                 nfreq++;
                 bcsav[nfreq] = omroot;
                 norm[nfreq] = dnor;
                 if (debug) printf("Accepted Root: %f Deter: %g %d\n",
                                    omroot/conv,dom1,nfreq);
                } 
              else
                {
                 if (debug) printf("Rejected Root: %f Deter: %g\n",
                                    omroot/conv,dom1);
                }
             }
           else
             {
              dnor *= pow( 1.0/fabs(func((om1+om2)/2.0)),1.0/nr); 
             }
           omega += del;
           om1 = om2;
           dom1 = dom2;
          }
     if (nfreq == 0) return;
     if (n == 0) 
       {
        effmass(bcsav,nfreq,spot,freq,weight);
       }
     else
       {
        for (nf = 1; nf <= nfreq; nf++)
           {
            freq[nf] = bcsav[nf] / conv;
            weight[nf] = 0.0;
           }
       }
     for (nf = 1; nf <=nfreq; nf++)
        {
         if (weight[nf] < 1.0e30)
           {
            reson[*NRF] = freq[nf];
            alpha[*NRF] = weight[nf];
            parity[*NRF] = iterm;
            order[*NRF] = nf;
            *NRF = *NRF + 1;
            printf(" %f %f %d %d %d \n",freq[nf],weight[nf],n,iterm,nf);
           }
        }
    }
}

double findroot(func,x1,x2,xacc)
double (*func)();
double x1,x2,xacc;
{
 int j;
 double dx,f,fmid,xmid,rtb;

 f = (*func)(x1);
 fmid = (*func)(x2);
 if (f * fmid >= 0.0) 
   {
    printf("Root must be bracketed for bisection\n");
    exit(1);
   }
 rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2);
 for (j=1;j<=100;j++) 
    {
     fmid = (*func)(xmid = rtb + (dx *= 0.5));
     if (fmid <= 0.0) rtb = xmid;
     if (fabs(dx) < xacc || fmid == 0.0) return rtb;
    }
 printf("Too many bisections\n");
 exit(1);
}

double func(omega)
double omega;
{
 int nd;
 double det;
 extern int nr;
 extern double d[41][41];
 void coefs(), deter();

 nd = 40;
 coefs(omega);
 deter(d,&det,nr,nd);
 return(det);
}


void deter(a,det,n,md)
double a[41][41],*det;
int n,md;
{
 static int i1,i3,i,j,j2;
 static double sum;

 if (n == 1) goto eighteen;
 *det = 1;
 i1 = 1;
 one:
 i3 = i1;
 sum = fabs(a[i1][i1]);
 for (i=i1; i<=n; i++)
    {
     if (sum - fabs(a[i][i1]) < 0) goto two; else goto three;
     two:
     i3 = i;
     sum = fabs(a[i][i1]);
     three:
     ;
    }
 if ((i3 - i1) == 0) goto six; else goto four;
 four:
 for (j=1; j<=n; j++)
    {
     sum = -a[i1][j];
     a[i1][j] = a[i3][j];
     a[i3][j] = sum;
    }
 six:
 i3 = i1 + 1;
 for (i=i3; i<=n; i++)
    a[i][i1] = a[i][i1] / a[i1][i1];
 j2 = i1 - 1;
 if (j2 == 0) goto eleven; else goto eight;
 eight:
 for (j=i3; j<=n; j++)
    for (i=1; i<=j2; i++)
       a[i1][j] = a[i1][j] - a[i1][i] * a[i][j];
 eleven:
 j2 = i1;
 i1++;
 for (i=i1; i<=n; i++)
    for (j=1; j<=j2; j++)
       a[i][i1] = a[i][i1] - a[i][j] * a[j][i1];
 if ((i1 - n) == 0) goto fourteen; else goto one;
 fourteen:
 i3 = 1;
 j2 = n / 2;
 if ((2*j2 - n) == 0) goto sixteen; else goto fifteen;
 fifteen:
 i3 = 0;
 *det = a[n][n];
 sixteen:
 for (i=1; i<=j2; i++)
    {
     j = n - i + i3;
     *det = *det * a[i][i] * a[j][j];
    }
 return;
 eighteen:
 *det = a[1][1];
}


void coefs(om)
double om;
{
 static int i,j,k,l;
 static double en,ens,enm,oms,r,r1,r2,t1,t2,t3,t4;
 extern int n,nz,iterm;
 extern double sd,v,h,aa[41],aas[41],bas[41],sba[41],ba[41],das[41],
        sda[41],da[41],bk[41],dk[41],bp[41],dp[41],bj[41],dj[41],ab2[41],
        bbs[41],sbb[41],bb[41],dbs[41],sdb[41],db[41],ss[41],cs[41],
        csd[41],b[41][41],bh[41][41],aba[41][41],cba[41][41],abas[41][41],
        bbas[41][41];
 extern double rj[51][6];
 void eval();
 double c(),s();

 sd = 1;
 en = n;
 ens = n * n;
 enm = ens - 1;
 oms = om * om;
 r = (1 - 2 * v) / (2 * (1 - v));
 for (i=1; i<=nz; i++)
    {
     if (iterm == 0) aa[i] = PI * (i - 1) / h;
     if (iterm == 1) aa[i] = PI * (2*i-1) / (2*h);
     aas[i] = aa[i] * aa[i];
     bas[i] = oms - aas[i];
     sba[i] = (bas[i] < 0 ? -1 : +1);
     ba[i] = sqrt(fabs(bas[i]));
     das[i] = oms * r - aas[i];
     sda[i] = (das[i] < 0 ? -1 : +1);
     da[i] = sqrt(fabs(das[i]));
     eval(n,ba[i],&bk[i],&bp[i],&bj[i],sba[i]);
     eval(n,da[i],&dk[i],&dp[i],&dj[i],sda[i]);
     if (n != 0)
       {
        a[i] = h * (2 * bp[i] * dj[i] * (bas[i] - aas[i]) * 
               (enm + (2 * aas[i] - oms) / 2 - dk[i]) + 4 * 
               aas[i] * dp[i] * bj[i] * (enm - bas[i] - bk[i]));
        cp[i] = 2 * h * en * bj[i] * (-2 * dp[i] * bk[i] + 
                dj[i] * (enm + 0.5 * (2 * aas[i] - oms) - dk[i]));
        ah[i] = h * (-2 * (bas[i] - aas[i]) * en * bp[i] * dj[i] * dk[i] - 
                4 * aas[i] * en * dp[i] * bj[i] * bk[i]);
        ch[i] = 2 * h  * bj[i] * (2 * dp[i] * (enm - bas[i] / 2 - bk[i]) -
                ens * dj[i] * dk[i]);
        sd *= (ch[i] < 0 ? -1 : +1);
       }
     else
       {
        a[i] = h * (2 * bp[i] * dj[i] * (bas[i] - aas[i]) * 
               (enm + (2 * aas[i] - oms) / 2 - dk[i]) + 4 * 
               aas[i] * das[i] * dp[i] * bj[i] * (-bp[i] / bj[i] - 1));
       }
    }
 for (j=1; j<=nr; j++)
    {
     ab[j] = rj[j][n+1];
     ab2[j] = ab[j] * ab[j];
     bbs[j] = oms - ab2[j];
     sbb[j] = (bbs[j] < 0 ? -1 : +1);
     bb[j] = sqrt(fabs(bbs[j]));
     dbs[j] = oms * r - ab2[j];
     sdb[j] = (dbs[j] < 0 ? -1 : +1);
     db[j] = sqrt(fabs(dbs[j]));
     if (iterm != 1)
       {
        ss[j]  = s(bb[j] * h, sbb[j]) * s(db[j] * h, sdb[j]) * dbs[j] * h * h;
        cs[j]  = c(db[j] * h, sdb[j]) * s(bb[j] * h, sbb[j]) * h;
        csd[j] = c(bb[j] * h, sbb[j]) * s(db[j] * h, sdb[j]) * h * dbs[j];
       }
     else
       {
        ss[j]  = -c(bb[j] * h, sbb[j]) * c(db[j] * h, sdb[j]);
        cs[j]  =  c(bb[j] * h, sbb[j]) * s(db[j] * h, sdb[j]) * h;
        csd[j] =  c(db[j] * h, sdb[j]) * s(bb[j] * h, sbb[j])* h * bbs[j];
       }
     if (n == 0) ao2[j] = 0.5;
     if (n > 0) ao2[j] = (ab2[j] - ens) / (2 * ab2[j]);
     bba[j] = ao2[j] * (2 * (ab2[j] - bbs[j]) * (ab2[j] - bbs[j]) * 
              cs[j] + 8 * ab2[j] * csd[j]);
    }
 for (i=1; i<=nz; i++)
    {
     for (j=1; j<=nr; j++)
        {
         t1 = 1 / (dbs[j] - aas[i]);
         sd *= (t1 < 0 ? -1 : +1);
         t2 = 1 / (bbs[j] - aas[i]);
         sd *= (t2 < 0 ? -1 : +1);
         t3 = 1 / (ab2[j] - bas[i]);
         t4 = 1 / (ab2[j] - das[i]);
         b[i][j] = (4 * (ab2[j] - bbs[j]) * (ens + (2 * dbs[j] - oms) / 2) 
                 * t1 + 8 * (ens - ab2[j]) * t2 * bbs[j]) * ss[j];
         if (n != 0)
           {
            bh[i][j] = -en * (-4 * (ab2[j] - bbs[j]) * t1 - 
                       8 * t2 * bbs[j]) * ss[j];
            aba[j][i] = 2 * (bas[i] - aas[i]) * bp[i] * (2 * das[i] - oms) * 
                        t4 * dp[i] + 8 * aas[i] * dp[i] * t3 * bas[i] * bp[i];
            cba[j][i] = 2 * en * bj[i] * (2 * das[i] - oms) * dp[i] * t4;
           }
         else
           {
            aba[j][i] = 2 * (bas[i] - aas[i]) * bp[i] * (2 * das[i] - oms) * 
                        t4 * dp[i] * das[i] + 8 * aas[i] * dp[i] * das[i] * 
                        t3 * bp[i] * bas[i];

           }
        }
    }
 if (iterm == 0) a[1] *= 2;
 if (n == 0)
   {
    for (i=1; i<=nz; i++)
       {
        as[i] = a[i];
        sd *= (as[i] < 0 ? -1 : +1);
        for (j=1; j<=nr; j++)
           {
            bs[i][j] = b[i][j];
            abas[j][i] = aba[j][i];
           }
       }
    for (j=1; j<=nr; j++)
       {
        bbas[j][j] = bba[j];
        for (k=1; k<=nr; k++)
           {
            if (k != j) bbas[j][k] = 0;
           }
       }
   }
 else
   {
    if (iterm == 0)
      {
       cp[1] *= 2;
       ah[1] *= 2;
       ch[1] *= 2;
      }
    for (i=1; i<=nz; i++)
       {
        r1 = cp[i] / ch[i];
        r2 = ah[i] / ch[i];
        as[i] = a[i] - ah[i] * r1;
        sd *= (as[i] < 0 ? -1 : +1);
        for (j=1; j<=nr; j++)
           {
            bs[i][j] = b[i][j] - bh[i][j] * r1;
            abas[j][i] = aba[j][i] - cba[j][i] * r2;
           }
       }
    for (j=1; j<=nr; j++)
       {
        bbas[j][j] = bba[j];
        for (k=1; k<=nr; k++)
           {
            if (k != j) bbas[j][k] = 0;
            for (l=1; l<=nz; l++)
               bbas[j][k] -= cba[j][l] * bh[l][k] / ch[l];
           }
       }
   }
 for (j=1; j<=nr; j++)
    {
     for (k=1; k<=nr; k++)
        {
         d[j][k] = bbas[j][k];
         for (l=1; l<=nz; l++)
            d[j][k] -= abas[j][l] * bs[l][k] / as[l];
        }
    }
 for (i=1; i<=nr; i++)
    for (j=1; j<=nr; j++)
        d[i][j] *= dnor;
}


void eval(n,x,fk,fp,fj,s)
int n;
double x,*fk,*fp,*fj,s;
{
 static int one;
 static double j[7],xton;
 double fn();
 void besj(),bessi();

 if (n == 0) goto five;
 if (x == 0) goto ten;
 if (s > 0) besj(&n,x,j);
 if (s < 0) bessi(n,x,j);
 xton = pow(x,n);
 *fk = x * j[n] / j[n+1] - n - 1;
 *fp = (x * j[n] - n * j[n+1]) / xton;
 *fj = j[n+1] / xton;
 return;
 five:
 one = 1;
 if (s > 0) besj(&one,x,j);
 one = 1;
 if (s < 0) bessi(one,x,j);
 *fj = j[1];
 if (x == 0) goto seven;
 *fk = -x * j[2] / j[1] * s - 1;
 *fp = -j[2] / x;
 return;
 seven:
 *fk = -1;
 *fp = -0.5;
 return;
 ten:
 *fk = n - 1;
 *fj = 1 / (fn(n) * pow(2,n));
 *fp = *fj * n;
}


void bessi(n,x,bj)
int n;
double x,bj[7];
{
 static int ier;
 static double bi;
 void besi();

 besi(x,n,&bi,&ier);
 bj[n+1] = bi;
 besi(x,n-1,&bi,&ier);
 bj[n] = bi;
}


double s(x,s1)
double x,s1;
{
 static double a,t;
 
 t = fabs(x);
 if ((s1 < 0) && (x < 30)) a = (1 - exp(-2*t)) / (2 * t);
 if ((s1 < 0) && (x >= 30)) a = 0.5 / t;
 if (s1 > 0) a = sin(x) / x;
 if (s1 == 0) a = 1;
 return (a);
}


double c(x,s1)
double x,s1;
{
 static double a,t;

 t = fabs(x);
 if ((s1 < 0) && (x < 30)) a = (1 + exp(-2*t)) / 2;
 if ((s1 < 0) && (x >= 30)) a = 0.5;
 if (s1 > 0) a = cos(x);
 if (s1 == 0) exit(17);
 return (a);
}


double fn(nn)
int nn;
{
 static int nf,i;
 nf = 1;
 for (i=2; i<=nn; i++) nf *= i;
 return((double)nf);
}


void besj(nmax,xp,j)
int *nmax;
double xp,*j;
{
 static int m,minus,ncount,nmp1,n,nu,np1,nr,limit,nlamb,lambda;
 static double ja[100],x,etol,d,eps,sum,d1,r,s,enmax,amax,en,alam;
 double tofy();

 ncount = 1;
 x = xp;
 etol = 1.0e-15;
 d = 11;
 if (x*(*nmax) < 0) minus = 1; else minus = 0;
 if (*nmax < 0) *nmax = - (*nmax);
 if (x < 0) x = -x;
 nmp1 = *nmax + 1;
 if (x < etol) goto five00;
 for (n=1; n<=nmp1; n++) ja[n] = 0;
 eps = 0.5e-11;
 sum = 1;
 d1 = 2.3026 * d + 1.3863;
 r = 0;
 enmax = *nmax;
 if (*nmax > 0) r = enmax * tofy(0.5 * d1 / enmax);
 s = 1.3591 * x * tofy(0.73576 * d1 / x);
 amax = s;
 if (r > s) amax = r;
 nu = 1.0 + amax;
 fifty:
 limit = nu / 2;
 n = 2 * limit;
 nlamb = 1;
 r = 0;
 s = 0;
 en = n;
 sixty:
 r = 1.0 / (2 * en / x - r);
 lambda = 1 + nlamb;
 nlamb = -nlamb;
 alam = lambda;
 s = r * (alam + s);
 nr = n + nmp1;
 if (n <= *nmax) ja[nr] = r;
 n--;
 en = n;
 if (n > 0) goto sixty;
 j[1] = sum / (1 + s);
 for (n=1; n <=(*nmax); n++)
    {
     np1 = n + 1;
     nr = n + nmp1;
     j[np1] = ja[nr] * j[n];
    }
 for (n=1; n<=nmp1; n++)
    {
     if (fabs((j[n]-ja[n])/j[n]) >= eps)
       {
        if (ncount <= 5)
          {
           ncount++;
           for (m=1; m<=nmp1; m++)
              ja[m] = j[m];
           nu += 5;
           goto fifty;
          }
       }
    }
 one00:
 if (!minus) goto two00;
 for (n=2; n<=nmp1; n+=2) j[n] = -j[n];
 two00:
 return;
 five00:
 j[1] = 1;
 if (*nmax == 0) goto two00;
 for (n=2; n<=nmp1; n++) j[n] = 0;
 goto two00;
}


double tofy(y)
double y;
{
 static double a,p,z;
 if (y <= 10)
   {
    p = (((0.000057941 * y - 0.00176148) * y + 0.0208645) * 
          y - 0.129013) * y + 0.85777;
    a = p * y + 1.0125;
   }
 else
   {
   z = log(y) - 0.775;
   p = 1.0/(1.0 + (0.775 - log(z)) / (1.0 + z));
   a = p * y / z;
   }
 return(a);
}


void besi(x,n,bi,ier)
int n,*ier;
double x,*bi;
{
 static int i,k;
 static double tol,xx,term,fi,fk,fn;

 *ier = 0;
 *bi = 1;
 if (n < 0) goto one50;
 if (n == 0) goto fifteen;
 if (n > 0) goto ten;
 ten:
 if (x < 0) goto one60;
 if (x >= 0) goto twenty;
 fifteen:
 if (x < 0) goto one60;
 if (x == 0) goto seventeen;
 if (x > 0) goto twenty;
 seventeen:
 return;
 twenty:
 tol = 1.0e-9;
 if ((x-12) <= 0) goto forty; else goto thirty;
 thirty:
 if ((x-n) <= 0) goto forty; else goto one10;
 forty:
 xx = 0.5 * x;
 fifty:
 term = 1;
 if (n <= 0) goto seventy; else goto fiftyfive;
 fiftyfive:
 for (i=1; i<=n; i++)
    {
     fi = i;
     if ((fabs(term) - 1.0e-36) < 0)
       {
        *ier = 3;
        *bi = 0;
        return;
       }
     else
       term *= xx / fi;
    }
 seventy:
 *bi = term;
 xx *= xx;
 for (k=1; k<=1000; k++)
    {
     if (fabs(term) - fabs(*bi * tol) > 0)
       {
        fk = k * (n + k);
        term *= xx / fk;
        *bi += term;
       }
     else
       goto one00;
    }
 one00:
 *bi *= exp(-x);
 return;
 one10:
 fn = 4 * n * n;
 xx = 0.125 / x;
 term = 1;
 *bi = 1;
 for (k=1; k<=30; k++)
    {
     if (fabs(term) - fabs(*bi * tol) > 0)
       {
        fk = (2*k-1);
        fk *= fk;
        term *= xx * (fk - fn) / k;
        *bi += term;
       }
     else 
       goto one40;
    }
 goto forty;
 one40:
 *bi /= sqrt(2 * PI * x);
 return;
 one50:
 *ier = 1;
 return;
 one60:
 *ier = 2;
}


void vecmat(nf)
int nf;
{
 extern int nr;
 extern double d[41][41], vecb[41][2001];
 double a[41][40], b[40], x[40];
 double w[40],v[40][40];
 double B[41],RHS[41];
 double magn,sum,wmax,wmin;
 int i,j,k,l,n,m;
 void svbksb(), svdcmp();

 m = nr; n = nr - 1;
 i=1;
    {
     for (j=1;j<=m;j++)
        {
         b[j] = -d[j][i];
         l = 1;
         for (k=1;k<=m;k++)
            {
             if (k != i)
               {
                a[j][l] = d[j][k];
                l++;
               }
            }
        }
     svdcmp(a,m,n,w,v);
     wmax = 0.0;
     for (j=1;j<=n;j++) if (w[j] > wmax) wmax = w[j];
     wmin = wmax * 1.0e-12;
     for (j=1;j<=n;j++) if (w[j] < wmin) w[j] = 0.0;
     svbksb(a,w,v,m,n,b,x);
     l = 1;
     for (j=1;j<=m;j++)
        {
         if (j == i)
           B[j] = 1.0;
         else
           {
            B[j] = x[l];
            l++;
           }
        }
     magn = 0.0;
     for (j=1; j<=nr; j++)
        {
         sum = 0.0;
         for (k=1; k<=nr; k++)
            {
             if (dbgall) printf("%.2e ",d[j][k]);
             sum += d[j][k] * B[k];
            }
         magn += sum * sum;
         if (dbgall) printf("   ");
         if (dbgall) printf("%.2e",B[j]);
         if (dbgall) printf("   ");
         if (dbgall) printf("%.2e",sum);
         if (dbgall) printf("\n");
        }
     if (dbgall) printf("\n");
     magn = sqrt(magn);
    }
 if (debug) printf("magn = %.4e\n",magn);
 for (i=1; i<=nr; i++) vecb[i][nf] = B[i];
}


void svbksb(u, w, v, m, n, b, x)
double u[41][40], w[40], v[40][40];
int m, n;
double b[40], x[40];
{
 int jj,j,i;
 double s,tmp[40];
 
 for (j=1;j<=n;j++)
    {
     s=0.0;
     if (w[j])
       {
        for (i=1;i<=m;i++) s += u[i][j]*b[i];
        s /= w[j];
       }
     tmp[j]=s;
    }
 for (j=1;j<=n;j++)
    {
     s=0.0;
     for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
     x[j]=s;
    }
}


static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
 
double pythag(a, b)
double a, b;
{
 double absa,absb;

 absa=fabs(a);
 absb=fabs(b);
 if (absa > absb)
   return absa*sqrt(1.0+SQR(absb/absa));
 else
   return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}


static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
        (dmaxarg1) : (dmaxarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
        (iminarg1) : (iminarg2))

void svdcmp(a, m, n, w, v)
double a[41][40];
int m, n; 
double w[40], v[40][40];
{
 double pythag();
 int flag,i,its,j,jj,k,l,nm;
 double anorm,c,f,g,h,s,scale,x,y,z,rv1[40];
 
 g=scale=anorm=0.0;
 for (i=1;i<=n;i++)
    {
     l=i+1;
     rv1[i]=scale*g;
     g=s=scale=0.0;
     if (i <= m)
       {
        for (k=i;k<=m;k++) scale += fabs(a[k][i]);
        if (scale)
          {
           for (k=i;k<=m;k++)
              {
               a[k][i] /= scale;
               s += a[k][i]*a[k][i];
              }
           f=a[i][i];
           g = -SIGN(sqrt(s),f);
           h=f*g-s;
           a[i][i]=f-g;
           for (j=l;j<=n;j++)
              {
               for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
               f=s/h;
               for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
              }
           for (k=i;k<=m;k++) a[k][i] *= scale;
          }
       }
     w[i]=scale *g;
     g=s=scale=0.0;
     if (i <= m && i != n)
       {
        for (k=l;k<=n;k++) scale += fabs(a[i][k]);
        if (scale)
          {
           for (k=l;k<=n;k++)
              {
               a[i][k] /= scale;
               s += a[i][k]*a[i][k];
              }
           f=a[i][l];
           g = -SIGN(sqrt(s),f);
           h=f*g-s;
           a[i][l]=f-g;
           for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
           for (j=l;j<=m;j++)
              {
               for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
               for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
              }
           for (k=l;k<=n;k++) a[i][k] *= scale;
          }
       }
     anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
    } 
 for (i=n;i>=1;i--)
    {
     if (i < n)
       {
        if (g)
          {
           for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
           for (j=l;j<=n;j++)
              {
               for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
               for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
              }
          }
        for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
       }
     v[i][i]=1.0;
     g=rv1[i];
     l=i;
    }
 for (i=IMIN(m,n);i>=1;i--)
    {
     l=i+1;
     g=w[i];
     for (j=l;j<=n;j++) a[i][j]=0.0;
     if (g)
       {
        g=1.0/g;
        for (j=l;j<=n;j++)
           {
            for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
            f=(s/a[i][i])*g;
            for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
           }
        for (j=i;j<=m;j++) a[j][i] *= g;
       }
     else for (j=i;j<=m;j++) a[j][i]=0.0;
     ++a[i][i];
    }
 for (k=n;k>=1;k--)
    {
     for (its=1;its<=30;its++)
        {
         flag=1;
         for (l=k;l>=1;l--)
            {
             nm=l-1;
             if ((double)(fabs(rv1[l])+anorm) == anorm)
               {
                flag=0;
                break;
               }
             if ((double)(fabs(w[nm])+anorm) == anorm) break;
            }
         if (flag)
           {
            c=0.0;
            s=1.0;
            for (i=l;i<=k;i++)
               {
                f=s*rv1[i];
                rv1[i]=c*rv1[i];
                if ((double)(fabs(f)+anorm) == anorm) break;
                g=w[i];
                h=pythag(f,g);
                w[i]=h;
                h=1.0/h;
                c=g*h;
                s = -f*h;
                for (j=1;j<=m;j++)
                   {
                    y=a[j][nm];
                    z=a[j][i];
                    a[j][nm]=y*c+z*s;
                    a[j][i]=z*c-y*s;
                   }
               }
           }
         z=w[k];
         if (l == k)
           {
            if (z < 0.0)
              {
               w[k] = -z;
               for (j=1;j<=n;j++) v[j][k] = -v[j][k];
              }
            break;
           }
         if (its == 60) printf("no convergence in 60 svdcmp iterations\n");
         x=w[l];
         nm=k-1;
         y=w[nm];
         g=rv1[nm];
         h=rv1[k];
         f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
         g=pythag(f,1.0);
         f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
         c=s=1.0;
         for (j=l;j<=nm;j++)
            {
             i=j+1;
             g=rv1[i];
             y=w[i];
             h=s*g;
             g=c*g;
             z=pythag(f,h);
             rv1[j]=z;
             c=f/z;
             s=h/z;
             f=x*c+g*s;
             g = g*c-x*s;
             h=y*s;
             y *= c;
             for (jj=1;jj<=n;jj++)
                {
                 x=v[jj][j];
                 z=v[jj][i];
                 v[jj][j]=x*c+z*s;
                 v[jj][i]=z*c-x*s;
                }
             z=pythag(f,h);
             w[j]=z;
             if (z)
               {
                z=1.0/z;
                c=f*z;
                s=h*z;
               }
             f=c*g+s*y;
             x=c*y-s*g;
             for (jj=1;jj<=m;jj++)
                {
                 y=a[jj][j];
                 z=a[jj][i];
                 a[jj][j]=y*c+z*s;
                 a[jj][i]=z*c-y*s;
                }
            }
         rv1[l]=0.0;
         rv1[k]=f;
         w[k]=x;
        }
    }
}


void vectors(nf)
int nf;
{
 static int i,j;
 static double btemp;
 extern int n,nz,nr;
 extern double bs[41][41],bh[41][41],as[41],ah[41],ch[41],
               veca[41][2001],vecb[41][2001],vecc[41][2001];

 for (i=1; i<=nz; i++)
    {
     btemp = 0;
     for (j=1; j<=nr; j++)
        {
         btemp += bs[i][j] * vecb[j][nf];
        }
     veca[i][nf] = -btemp / as[i];
    }
 if (n != 0)
   {
    for (i=1; i<=nz; i++)
       {
        btemp = 0;
        for (j=1; j<=nr; j++)
           {
            btemp += bh[i][j] * vecb[j][nf];
           }
        vecc[i][nf] = -(ah[i] * veca[i][nf] + btemp) / ch[i];
       }
   }
}


void vals()
{
 static int i;
 static double bbh,dbh,sab;
 extern int n,nz,nr;
 extern double h,ba[41],da[41],ab[41],sba[41],sda[41],db[41],sdb[41],sbb[41];
 extern double psiba[41],qjbba[41],qjpba[41],psida[41],qjpda[41],
               qjbda[41],psiab[41],qjpab[41],qjbab[41],cosdbh[41],
               cosbbh[41],cosaah[41],sndbh[41],snbbh[41],snaah[41];
 void eval();
 double c();

 for (i=1;i<=nz; i++)
    {
     eval(n,ba[i],&psiba[i],&qjpba[i],&qjbba[i],sba[i]);
     eval(n,da[i],&psida[i],&qjpda[i],&qjbda[i],sda[i]);
    }
 sab = 2;
 for (i=1;i<=nr; i++)
    {
     eval(n,ab[i],&psiab[i],&qjpab[i],&qjbab[i],sab);
    }
 for (i=1;i<=nr; i++)
    {
     dbh = db[i] * h;
     cosdbh[i] = c(dbh,sdb[i]);
     bbh = bb[i] * h;
     cosbbh[i] = c(bbh,sbb[i]);
     sndbh[i] = s(dbh,sdb[i]);
     snbbh[i] = s(bbh,sbb[i]);
    }
  for (i=1;i<=nz; i++)
    {
     if (aa[i] != 0) snaah[i] = sin(aa[i] * h) / (aa[i] * h);
     if (aa[i] == 0) snaah[i] = 1;
     cosaah[i] = cos(aa[i] * h);
    }
}


void valsrz(r,z)
double r,z;
{
 static int i;
 static double sab,bar,dar,abr,t,ed,eb,dbz,bbz,q;
 extern int n,nz,nr;
 extern double ba[41],sba[41],da[41],sda[41],ab[41];
 extern double psibar[41],qjbbar[41],qjpbar[41],psidar[41],qjpdar[41],
               qjbdar[41],psiabr[41],qjpabr[41],qjbabr[41],cosdbz[41],
               cosbbz[41],cosaaz[41],sndbz[41],snbbz[41],snaaz[41];
 void eval();
 double c(),s();

 for (i=1;i<=nz; i++)
    {
     bar = r * ba[i];
     eval(n,bar,&psibar[i],&qjpbar[i],&qjbbar[i],sba[i]);
     if (sba[i] < 0) qjpbar[i] *= exp(ba[i] * (r-1));
     if (sba[i] < 0) qjbbar[i] *= exp(ba[i] * (r-1));
     dar = r * da[i];
     eval(n,dar,&psidar[i],&qjpdar[i],&qjbdar[i],sda[i]);
     if (sda[i] < 0) qjpdar[i] *= exp(da[i] * (r-1));
     if (sda[i] < 0) qjbdar[i] *= exp(da[i] * (r-1));
    }
 sab = 2;
 for (i=1;i<=nr; i++)
    {
     abr = ab[i] * r;
     eval(n,abr,&psiabr[i],&qjpabr[i],&qjbabr[i],sab);
    }
 for (i=1;i<=nr; i++)
    {
     t = fabs(z);
     ed = exp(db[i] * (t - h));
     eb = exp(bb[i] * (t - h));
     dbz = db[i] * z;
     cosdbz[i] = c(dbz,sdb[i]);
     if (sdb[i] < 0) cosdbz[i] = cosdbz[i] * ed;
     bbz = bb[i] * z;
     cosbbz[i] = c(bbz,sbb[i]);
     if (sbb[i] < 0) cosbbz[i] = cosbbz[i] * eb;
     sndbz[i] = s(dbz,sdb[i]);
     if (sdb[i] < 0) sndbz[i] = sndbz[i] * ed;
     snbbz[i] = s(bbz,sbb[i]);
     if (sbb[i] < 0) snbbz[i] = snbbz[i] * eb;
    }
 for (i=1;i<=nz; i++)
    {
     q = aa[i] * z;
     if (q != 0) snaaz[i] = sin(q) / q;
     if (q == 0) snaaz[i] = 1;
     cosaaz[i] = cos(q);
    }
}


void displ0(nf,xr,xz,r,z)
int nf;
double *xr,*xz,r,z;
{
 static int i,j;
 static double u5,w4,t1,u1,u2,w2,u6,w1,u3,w3,u4;
 extern int nz,nr;
 extern double aas[41],bas[41],qjpda[41],qjpbar[41],qjpba[41],qjpdar[41],
               snaaz[41],snaah[41],cosaaz[41],cosaah[41],
               veca[41][2001],vecb[41][2001],vecc[41][2001],das[41],
               qjbbar[41],qjbdar[41],bbs[41],ab2[41],qjpabr[41],
               qjbab[41],cosdbh[41],snbbz[41],cosbbh[41],sndbz[41],
               dbs[41],sndbh[41],coshhz[41],snbbh[41],cosdbz[41],
               qjbabr[41],cosbbz[41];

 u5 = 0;
 w4 = 0;
 for (i=1;i<=nz; i++)
    {
     t1 = aas[i] - bas[i];
     u1 = 2 * aas[i] * qjpda[i] * qjpbar[i] - t1 * qjpba[i] * qjpdar[i];
     if (iterm == 1) u4 = z / h * snaaz[i] / snaah[i];
     if (iterm == 0) u4 = cosaaz[i] / cosaah[i];
     u6 = r * das[i] * u1 * u4;
     u5 = u6 * veca[i][nf] + u5;
     w1 = 2 * das[i] * qjpda[i] * qjbbar[i] + t1 * qjpba[i] * qjbdar[i];
     if (iterm == 1) w3 = -cosaaz[i] / snaah[i] / h;
     if (iterm == 0) w3 = aas[i] * z * snaaz[i] / cosaah[i];
     w4 = w1 * w3 * veca[i][nf] + w4;
    }
 u3 = 0;
 w3 = 0;
 for (j=1;j<=nr; j++)
    {
     t1 = bbs[j] - ab2[j];
     u1 = ab2[j] * r * qjpabr[j] / qjbab[j] ;
     if (iterm == 1) u2 = 2 * bbs[j] * z * cosdbh[j] * snbbz[j] - 
                          t1 * z * cosbbh[j] * sndbz[j];
      if (iterm == 0) u2 = 2 * dbs[j] * h * sndbh[j] * cosbbz[j] - 
                           t1 * h * snbbh[j] * cosdbz[j];
      u3 = u1 * u2 * vecb[j][nf] + u3;
      w1 = qjbabr[j] / qjbab[j];
      if (iterm == 1) w2 = -2 * ab2[j] * cosdbh[j] * cosbbz[j] -
                           t1 * cosbbh[j] * cosdbz[j];
      if (iterm == 0) w2 = (2 * ab2[j] * sndbh[j] * snbbz[j] +
                           t1 * snbbh[j] * sndbz[j]) * dbs[j] * z * h;
      w3 = w1 * w2 * vecb[j][nf] + w3;
    }
 *xr = u5 + u3;
 *xz = w4 + w3;
}


void energy(om,nf,dxz1)
int nf;
double om,*dxz1;
{
 static int i,j;
 static double t,dz,dz2,dr,dr2,r1,r2,r3,omega,ener,dxz,z,r,xr,om2rad2;
 static double denstmass,xz0,dener,ratio,ratsq,xz[51][51],dvol[51][51];
 extern int ndz,ndr;
 extern double h,rad,dens,shear,speed,tmass,veca[41][2001],vecb[41][2001];
 void coefs(),vals(),valsrz(),displ0();

 coefs(om);
 vecmat(nf);
 vectors(nf);
 vals();
 t = ndz;
 dz = 2 * h / t;
 dz2 = dz / 2;
 dr = ndr;
 dr = 1 / dr;
 dr2 = dr / 2;
 r3 = rad * rad * rad;
 omega = om / rad * speed;
 om2rad2 = omega * rad;
 om2rad2 *= om2rad2;
 denstmass = dens / tmass;
 z = -h - dz2;
 ener = 0;
 dxz = 0;
 for (i=1; i<=ndz; i++)
    {
     z += dz;
     r = -dr2;
     for (j=1; j<=ndr; j++)
        {
         r += dr;
         valsrz(r,z);
         displ0(nf,&xr,&xz[j][i],r,z);
         r1 = r - dr2;
         r2 = r + dr2;
         dvol[j][i] = dz * PI * (r2 * r2 - r1 * r1) * r3;
         ener += 0.5 * dvol[j][i] * dens * xr * xr * om2rad2; 
         dxz += xz[j][i] * dvol[j][i] * denstmass;
        }
    }
 for (i=1; i<=ndz; i++)
    {
     for (j=1; j<=ndr; j++)
        {
         xz0 = xz[j][i] - dxz;
         dener = 0.5*dvol[j][i] * dens * xz0 * xz0 * om2rad2;
         ener += dener;
        }
    }
 ratio = ener;
 ratsq = sqrt(ratio);
 *dxz1 = dxz / ratsq;
 for (i=1; i<=nz; i++) veca[i][nf] /= ratsq;
 for (j=1; j<=nr; j++) vecb[j][nf] /= ratsq;
}


void effmass(bcsav,nfreq,spot,freq,weight)
int nfreq;
double *bcsav,*freq,*weight,spot;
{
 static int i,nf;
 static double tmrad2,w2z,z,dxz1,omega,dr,r,zmove,xr,xz,area,qintens,xz1,ener;
 extern double h,rad,shear,speed,dens,tmass,dnor,norm[2001];
 void energy();

 tmass = dens * 2 * h * rad * rad * rad * PI;
 tmrad2 = tmass * rad * rad;
 w2z = (spot*spot) / (2*rad*rad);
 z = h;
 for (nf=1; nf<=nfreq; nf++)
    {
     dnor = norm[nf];
     energy(bcsav[nf],nf,&dxz1);
     omega = bcsav[nf] / rad * speed;
     dr = 1.0e-3;
     r = -dr / 2;
     zmove = 0;
     for (i=1; i<=1000; i++)
        {
         r += dr;
         valsrz(r,z);
         displ0(nf,&xr,&xz,r,z);
         area = PI * (SQR(r+dr/2) - SQR(r-dr/2));
         qintens = exp(-r * r / w2z) / PI / w2z;
         xz1 = xz - dxz1;
         zmove += xz1 * qintens * area;
        }
     ener = 0.5 * tmrad2 * zmove * zmove * omega * omega;
     weight[nf] = 1.0 / ener;
     freq[nf] = omega / (2 * PI);
    }
}
void opening_message(char name[]){
    printf("\n\n\n");
    printf("--------------------------------------------------------------\n");
    printf("\n");
    printf("                        %s\n",name);
    printf("      a program to calculate resonant frequencies\n");
    printf("                            &\n");
    printf("             corresponding effective masses\n");
    printf("                            of\n");
    printf("             right circular solid cylinders\n");
    printf("                                             LIGO project 1997\n");
    printf("--------------------------------------------------------------\n");
    printf("USAGE: %s \n", name);
    printf("   OR: %s -file input_file_name \n", name);
    printf("--------------------------------------------------------------\n");
    printf("\n\n");
}