#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	CPP/
#	CPP/Dp/
#	CPP/Dp/Drivers/
#	CPP/Dp/Drivers/Makefile
#	CPP/Dp/Drivers/data1
#	CPP/Dp/Drivers/data2
#	CPP/Dp/Drivers/driver.cpp
#	CPP/Dp/Drivers/res1
#	CPP/Dp/Drivers/res2
#	CPP/Dp/Src/
#	CPP/Dp/Src/bsslr.cpp
#	CPP/Dp/Src/bsslr.h
#	CPP/Dp/Src/mathur.cpp
#	CPP/Dp/Src/mathur.h
#	CPP/Dp/Src/mbsslr.cpp
#	CPP/Dp/Src/mbsslr.h
#	CPP/Dp/Src/mcnr.cpp
#	CPP/Dp/Src/mcnr.h
#	CPP/Dp/Src/mmathur.cpp
#	CPP/Dp/Src/mmathur.h
#	Doc/
#	Doc/Read.me
# This archive created: Tue Jul 23 15:16:31 2002
export PATH; PATH=/bin:$PATH
if test ! -d 'CPP'
then
	mkdir 'CPP'
fi
cd 'CPP'
if test ! -d 'Dp'
then
	mkdir 'Dp'
fi
cd 'Dp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
include makefile.inc

OBJS = driver.o bsslr.o mathur.o mbsslr.o mcnr.o mmathur.o

all: res1 res2

res1: $(OBJS) data1
	$(CPP) $(CPPOPTS) -o driver $(OBJS)
	driver < data1 > res1

res2: $(OBJS) data2
	$(CPP) $(CPPOPTS) -o driver $(OBJS)
	driver < data2 > res2

driver.o: driver.cpp $(INCSRC)/bsslr.cpp $(INCSRC)/mathur.cpp \
          $(INCSRC)/mbsslr.cpp $(INCSRC)/mcnr.cpp $(INCSRC)/mmathur.cpp
	$(CPP) $(CPPOPTS) -c -I$(INCSRC) driver.cpp

bsslr.o: $(INCSRC)/bsslr.cpp
	$(CPP) $(CPPOPTS) -c $(INCSRC)/bsslr.cpp

mathur.o: $(INCSRC)/mathur.cpp
	$(CPP) $(CPPOPTS) -c $(INCSRC)/mathur.cpp

mbsslr.o: $(INCSRC)/mbsslr.cpp
	$(CPP) $(CPPOPTS) -c $(INCSRC)/mbsslr.cpp

mcnr.o: $(INCSRC)/mcnr.cpp
	$(CPP) $(CPPOPTS) -c $(INCSRC)/mcnr.cpp

mmathur.o: $(INCSRC)/mmathur.cpp
	$(CPP) $(CPPOPTS) -c $(INCSRC)/mmathur.cpp
SHAR_EOF
fi # end of overwriting check
if test -f 'data1'
then
	echo shar: will not over-write existing file "'data1'"
else
cat << "SHAR_EOF" > 'data1'
e
3
c
f
s
SHAR_EOF
fi # end of overwriting check
if test -f 'data2'
then
	echo shar: will not over-write existing file "'data2'"
else
cat << "SHAR_EOF" > 'data2'
o
3
r
d
m
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.cpp'
then
	echo shar: will not over-write existing file "'driver.cpp'"
else
cat << "SHAR_EOF" > 'driver.cpp'
/*            Standard include libraries     */
#include   <cstdio>
#include   <cstdlib>
#include   <iostream>
#include   <fstream>
#include   <cmath>

/*       Include libraries  part of Mathieu functions Library  */
#include   "mathur.h"
#include   "bsslr.h"
#include   "mmathur.h"
#include   "mbsslr.h"
/*
  The Author would be happy to help in any way he can
  in the use of these routins.

  Contact Address:
  CERI, KACST, P.O. Box 6086,
  Riyadh 11442,
  Saudi Arabia,
  Fax:966+1+4813764
  Email: alhargan@kacst.edu.sa

  Legal Matters:
  The program is here as is, no guarantees are given, stated or implied,
  you may use it at your own risk.

  Last modified 20/5/1999

*/

struct point {
   double x;
   double y;
};
/*
  This a driver for the Mathieu functions library
  The following files must be compiled and linked:
  MATHUR.CPP     Mathieu functions
  MMATHUR.CPP    Modified Mathieu functions
  MCNR.CPP       Mathieu charactristic Numbers
  BSSLR.CPP      Bessel functions
  MBSSLR.CPP     Modified Bessel functions
*/
    double Radial(char typ,int n,double h,double u,int kind, double mdf);
    double Circumf(char typ,int n,double h,double v,int kind, double mdf);

/*----------------- Main Porgram --------------------------*/
int  main (void)
{
  point sn1[6];
  int i;
  int n,M;
  double h,v;
  double t,tm;
  double pi;
  char typ,r,mdc,kd;
  short kind=1,md=0;
  char num[]="0123456789";
  char ordc[]="00";
  int n1,n2;
  pi=4*atan(1.0);

  std::cout <<"\n typ=";
  std::cin >>typ;
  std::cout <<" order n=";
  std::cin >>n;
  std::cout <<" Circumferential(c) or Radial(r): ";
  std::cin >>r;
  std::cout <<" kind: First(f) Second(d): ";
  std::cin >>kd;
  std::cout <<" Standard(s) or Modified(m): ";
  std::cin >>mdc;

  if(n<10) ordc[1]=num[n];
	   else
	   {
	    n1=n/10;
	    n2=n-n1*10;
	    ordc[0]=num[n1];
	    ordc[1]=num[n2];
	   }

  if (kd=='d') kind=2; else kind=1;
  if (mdc=='m') md=1; else md=0;

   M=5;

 printf("\n\n%4d \n",n);
 if(r=='r')
 {
  M=4;
  tm=0.0125;
  if(kind==2) tm=10*0.0125;
  for(t=tm;t<=2.5001+tm;t +=0.0125)
  {
    h=t*n;
    for(i=1;i<=M;i++) sn1[i].x=t;

    sn1[1].y=Radial(typ,n,h,0.25,kind,md);  // u=2.5
    sn1[2].y=Radial(typ,n,h,0.5,kind,md);
    sn1[3].y=Radial(typ,n,h,1,kind,md);
    sn1[4].y=Radial(typ,n,h,2,kind,md);

    printf("\n%3.3f   ",sn1[1].x);
    for (i=1;i<=M;i++)
     {
       printf("%+1.6f    ",sn1[i].y);
     }

     }
  }
 else
 {
  M=4;
  for(v=0;v<=2*pi;v +=0.0314)
  {
    for(i=1;i<=M;i++) sn1[i].x=v*180/pi;

    sn1[1].y=Circumf(typ,n,0.5*n,v,kind,md);  // h=0.5*n
    sn1[2].y=Circumf(typ,n,1.0*n,v,kind,md);
    sn1[3].y=Circumf(typ,n,2.0*n,v,kind,md);
    sn1[4].y=Circumf(typ,n,3.0*n,v,kind,md);

    printf("\n%3.3f   ",sn1[1].x);
    for (i=1;i<=M;i++)
     {
       printf("%+1.6f     ",sn1[i].y);
     }

     }

 }


  printf("\n");
  exit(0);
}  // End Main

double Radial(char typ,int n,double h,double u,int kind, double mdf)
{

  if(mdf==1) return MathuMZn(typ,n,h,u,kind);
  return MathuZn(typ,n,h,u,kind);
}
double Circumf(char typ,int n,double h,double v,int kind, double mdf)
{

  if(mdf==1) return MathuQn(typ,n,h,v,kind);
  return MathuSn(typ,n,h,v,kind);
}
SHAR_EOF
fi # end of overwriting check
if test -f 'res1'
then
	echo shar: will not over-write existing file "'res1'"
else
cat << "SHAR_EOF" > 'res1'

 typ= order n= Circumferential(c) or Radial(r):  kind: First(f) Second(d):  Standard(s) or Modified(m): 

   3 

0.000   +1.000000     +1.000000     +1.000000     +1.000000     
1.799   +0.996109     +0.997544     +1.001524     +1.012443     
3.598   +0.984464     +0.990178     +1.006065     +1.050001     
5.397   +0.965149     +0.977913     +1.013532     +1.113356     
7.196   +0.938305     +0.960768     +1.023769     +1.203622     
8.995   +0.904126     +0.938769     +1.036555     +1.322318     
10.795   +0.862863     +0.911956     +1.051601     +1.471312     
12.594   +0.814818     +0.880379     +1.068541     +1.652759     
14.393   +0.760346     +0.844107     +1.086933     +1.869006     
16.192   +0.699854     +0.803227     +1.106252     +2.122464     
17.991   +0.633795     +0.757848     +1.125888     +2.415457     
19.790   +0.562670     +0.708106     +1.145142     +2.750024     
21.589   +0.487025     +0.654169     +1.163224     +3.127680     
23.388   +0.407444     +0.596234     +1.179262     +3.549143     
25.187   +0.324551     +0.534539     +1.192297     +4.014008     
26.986   +0.239002     +0.469363     +1.201303     +4.520399     
28.785   +0.151486     +0.401027     +1.205193     +5.064589     
30.584   +0.062714     +0.329900     +1.202841     +5.640621     
32.384   -0.026582     +0.256398     +1.193111     +6.239946     
34.183   -0.115652     +0.180987     +1.174879     +6.851109     
35.982   -0.203741     +0.104185     +1.147079     +7.459526     
37.781   -0.290088     +0.026555     +1.108737     +8.047392     
39.580   -0.373935     -0.051289     +1.059026     +8.593758     
41.379   -0.454534     -0.128694     +0.997308     +9.074848     
43.178   -0.531156     -0.204965     +0.923192     +9.464617     
44.977   -0.603093     -0.279375     +0.836579     +9.735624     
46.776   -0.669669     -0.351172     +0.737717     +9.860193     
48.575   -0.730245     -0.419589     +0.627239     +9.811884     
50.374   -0.784228     -0.483852     +0.506195     +9.567210     
52.174   -0.831077     -0.543191     +0.376074     +9.107551     
53.973   -0.870310     -0.596857     +0.238809     +8.421134     
55.772   -0.901510     -0.644130     +0.096763     +7.504970     
57.571   -0.924329     -0.684336     -0.047308     +6.366546     
59.370   -0.938496     -0.716860     -0.190304     +5.025133     
61.169   -0.943821     -0.741158     -0.328866     +3.512486     
62.968   -0.940197     -0.756777     -0.459470     +1.872803     
64.767   -0.927604     -0.763358     -0.578543     +0.161796     
66.566   -0.906112     -0.760655     -0.682598     -1.555186     
68.365   -0.875880     -0.748544     -0.768370     -3.205985     
70.164   -0.837156     -0.727027     -0.832966     -4.715307     
71.963   -0.790277     -0.696241     -0.873995     -6.009072     
73.763   -0.735666     -0.656459     -0.889707     -7.019154     
75.562   -0.673829     -0.608094     -0.879094     -7.688185     
77.361   -0.605349     -0.551691     -0.841975     -7.974016     
79.160   -0.530883     -0.487927     -0.779046     -7.853465     
80.959   -0.451151     -0.417598     -0.691895     -7.324965     
82.758   -0.366934     -0.341611     -0.582970     -6.409828     
84.557   -0.279063     -0.260969     -0.455519     -5.151924     
86.356   -0.188407     -0.176755     -0.313483     -3.615730     
88.155   -0.095870     -0.090115     -0.161354     -1.882816     
89.954   -0.002375     -0.002234     -0.004014     -0.047020     
91.753   +0.091143     +0.085677     +0.153457     +1.791349     
93.553   +0.183751     +0.172409     +0.305970     +3.531820     
95.352   +0.274524     +0.256774     +0.448630     +5.080081     
97.151   +0.362558     +0.337624     +0.576921     +6.353814     
98.950   +0.446980     +0.413871     +0.686870     +7.287577     
100.749   +0.526959     +0.484510     +0.775192     +7.836395     
102.548   +0.601710     +0.548628     +0.839395     +7.977790     
104.347   +0.670510     +0.605421     +0.877847     +7.712187     
106.146   +0.732699     +0.654208     +0.889808     +7.061724     
107.945   +0.787689     +0.694436     +0.875414     +6.067670     
109.744   +0.834971     +0.725687     +0.835632     +4.786725     
111.543   +0.874118     +0.747678     +0.772179     +3.286583     
113.343   +0.904787     +0.760268     +0.687415     +1.641139     
115.142   +0.926727     +0.763446     +0.584211     -0.074266     
116.941   +0.939774     +0.757331     +0.465817     -1.787220     
118.740   +0.943853     +0.742163     +0.335712     -3.431953     
120.539   +0.938980     +0.718294     +0.197468     -4.952209     
122.338   +0.925257     +0.686176     +0.054614     -6.303179     
124.137   +0.902872     +0.646345     -0.089479     -7.452465     
125.936   +0.872089     +0.599415     -0.231698     -8.380179     
127.735   +0.833254     +0.546056     -0.369266     -9.078258     
129.534   +0.786780     +0.486986     -0.499799     -9.549202     
131.333   +0.733146     +0.422955     -0.621343     -9.804381     
133.132   +0.672891     +0.354730     -0.732386     -9.862118     
134.932   +0.606605     +0.283086     -0.831856     -9.745698     
136.731   +0.534924     +0.208790     -0.919099     -9.481457     
138.530   +0.458523     +0.132596     -0.993850     -9.097043     
140.329   +0.378107     +0.055232     -1.056189     -8.619943     
142.128   +0.294407     -0.022607     -1.106496     -8.076293     
143.927   +0.208169     -0.100262     -1.145395     -7.490003     
145.726   +0.120150     -0.177121     -1.173708     -6.882172     
147.525   +0.031111     -0.252614     -1.192400     -6.270769     
149.324   -0.058191     -0.326225     -1.202535     -5.670543     
151.123   -0.147007     -0.397483     -1.205233     -5.093100     
152.922   -0.234605     -0.465970     -1.201633     -4.547131     
154.722   -0.320270     -0.531315     -1.192860     -4.038717     
156.521   -0.403314     -0.593193     -1.180003     -3.571686     
158.320   -0.483078     -0.651324     -1.164095     -3.148001     
160.119   -0.558938     -0.705470     -1.146094     -2.768132     
161.918   -0.630305     -0.755429     -1.126880     -2.431411     
163.717   -0.696633     -0.801032     -1.107246     -2.136355     
165.516   -0.757419     -0.842143     -1.087895     -1.880946     
167.315   -0.812206     -0.878652     -1.069443     -1.662870     
169.114   -0.860586     -0.910468     -1.052418     -1.479715     
170.913   -0.902202     -0.937525     -1.037267     -1.329128     
172.712   -0.936746     -0.959769     -1.024358     -1.208944     
174.511   -0.963968     -0.977161     -1.013986     -1.117276     
176.311   -0.983668     -0.989674     -1.006374     -1.052586     
178.110   -0.995705     -0.997288     -1.001682     -1.013740     
179.909   -0.999990     -0.999994     -1.000004     -1.000032     
181.708   -0.996493     -0.997786     -1.001373     -1.011211     
183.507   -0.985240     -0.990670     -1.005763     -1.047483     
185.306   -0.966312     -0.978653     -1.013085     -1.109506     
187.105   -0.939845     -0.961755     -1.023186     -1.198374     
188.904   -0.906033     -0.940001     -1.035850     -1.315585     
190.703   -0.865122     -0.913431     -1.050789     -1.462992     
192.502   -0.817413     -0.882094     -1.067642     -1.642739     
194.301   -0.763257     -0.846058     -1.085973     -1.857162     
196.101   -0.703060     -0.805409     -1.105259     -2.108675     
197.900   -0.637272     -0.760256     -1.124895     -2.399611     
199.699   -0.566391     -0.710732     -1.144186     -2.732027     
201.498   -0.490961     -0.657002     -1.162349     -3.107472     
203.297   -0.411565     -0.599265     -1.178512     -3.526712     
205.096   -0.328825     -0.537755     -1.191724     -3.989407     
206.895   -0.243395     -0.472748     -1.200960     -4.493764     
208.694   -0.155961     -0.404564     -1.205136     -5.036160     
210.493   -0.067235     -0.333568     -1.203129     -5.610759     
212.292   +0.022051     -0.260176     -1.193800     -6.209153     
214.091   +0.111152     -0.184850     -1.176026     -6.820040     
215.890   +0.199309     -0.108105     -1.148735     -7.428998     
217.690   +0.285762     -0.030503     -1.110950     -8.018384     
219.489   +0.369754     +0.047346     -1.061832     -8.567406     
221.288   +0.450536     +0.124789     -1.000735     -9.052417     
223.087   +0.527377     +0.201135     -0.927252     -9.447472     
224.886   +0.599568     +0.275657     -0.841271     -9.725172     
226.685   +0.666431     +0.347605     -0.743018     -9.857823     
228.484   +0.727326     +0.416212     -0.633108     -9.818880     
230.283   +0.781657     +0.480704     -0.512567     -9.584664     
232.082   +0.828880     +0.540312     -0.382864     -9.136258     
233.881   +0.868511     +0.594283     -0.245908     -8.461497     
235.680   +0.900126     +0.641897     -0.104041     -7.556899     
237.480   +0.923378     +0.682477     +0.039999     -6.429389     
239.279   +0.937989     +0.715404     +0.183129     -5.097613     
241.078   +0.943765     +0.740132     +0.322000     -3.592690     
242.877   +0.940597     +0.756199     +0.453093     -1.958201     
244.676   +0.928458     +0.763246     +0.572836     -0.249310     
246.475   +0.907414     +0.761018     +0.677734     +1.469064     
248.274   +0.877620     +0.749385     +0.764507     +3.125022     
250.073   +0.839319     +0.728343     +0.830239     +4.643333     
251.872   +0.792845     +0.698021     +0.872512     +5.949741     
253.671   +0.738615     +0.658687     +0.889539     +6.975704     
255.470   +0.677131     +0.610745     +0.880273     +7.663193     
257.270   +0.608973     +0.554735     +0.844489     +7.969192     
259.069   +0.534793     +0.491326     +0.782838     +7.869481     
260.868   +0.455310     +0.421309     +0.696863     +7.361354     
262.667   +0.371301     +0.345586     +0.588971     +6.464956     
264.466   +0.283594     +0.265155     +0.462371     +5.223048     
266.265   +0.193059     +0.181095     +0.320970     +3.699131     
268.064   +0.100595     +0.094549     +0.169238     +1.974017     
269.863   +0.007126     +0.006703     +0.012040     +0.141054     
271.662   -0.086414     -0.081236     -0.145547     -1.699629     
273.461   -0.179091     -0.168057     -0.298431     -3.447413     
275.260   -0.269978     -0.252570     -0.441704     -5.007528     
277.059   -0.358173     -0.333625     -0.570824     -6.296920     
278.859   -0.442799     -0.410131     -0.681789     -7.249194     
280.658   -0.523022     -0.481077     -0.771275     -7.818272     
282.457   -0.598056     -0.545546     -0.836748     -7.980513     
284.256   -0.667175     -0.602728     -0.876532     -7.735195     
286.055   -0.729713     -0.651936     -0.889840     -7.103407     
287.854   -0.785082     -0.692609     -0.876766     -6.125528     
289.653   -0.832766     -0.724323     -0.838237     -4.857578     
291.452   -0.872334     -0.746789     -0.775933     -3.366808     
293.251   -0.903440     -0.759856     -0.692184     -1.726913     
295.050   -0.925828     -0.763510     -0.589840     -0.013271     
296.849   -0.939328     -0.757861     -0.472134     +1.701460     
298.649   -0.943862     -0.743145     -0.342538     +3.351096     
300.448   -0.939441     -0.719708     -0.204620     +4.878848     
302.247   -0.926164     -0.687995     -0.061917     +6.239290     
304.046   -0.904212     -0.648542     +0.082191     +7.399389     
305.845   -0.873848     -0.601956     +0.224575     +8.338630     
307.644   -0.835412     -0.548906     +0.362439     +9.048379     
309.443   -0.789314     -0.490108     +0.493380     +9.530637     
311.242   -0.736031     -0.426310     +0.615420     +9.796370     
313.041   -0.676097     -0.358280     +0.727025     +9.863594     
314.840   -0.610102     -0.286790     +0.827101     +9.755392     
316.639   -0.538679     -0.212610     +0.914974     +9.497987     
318.438   -0.462501     -0.136495     +0.990359     +9.119000     
320.238   -0.382271     -0.059173     +1.053321     +8.645957     
322.037   -0.298720     +0.018658     +1.104225     +8.105085     
323.836   -0.212592     +0.096338     +1.143685     +7.520424     
325.635   -0.124646     +0.173251     +1.172512     +6.913226     
327.434   -0.035640     +0.248826     +1.191667     +6.301623     
329.233   +0.053667     +0.322544     +1.202210     +5.700524     
331.032   +0.142526     +0.393932     +1.205258     +5.121693     
332.831   +0.230202     +0.462569     +1.201949     +4.573961     
334.630   +0.315982     +0.528081     +1.193412     +4.063532     
336.429   +0.399175     +0.590142     +1.180737     +3.594341     
338.228   +0.479121     +0.648470     +1.164960     +3.168435     
340.028   +0.555193     +0.702823     +1.147043     +2.786351     
341.827   +0.626802     +0.752998     +1.127871     +2.447471     
343.626   +0.693398     +0.798826     +1.108241     +2.150347     
345.425   +0.754477     +0.840168     +1.088860     +1.892981     
347.224   +0.809578     +0.876912     +1.070348     +1.673070     
349.023   +0.858293     +0.908969     +1.053240     +1.488201     
350.822   +0.900259     +0.936268     +1.037984     +1.336017     
352.621   +0.935168     +0.958757     +1.024954     +1.214339     
354.420   +0.962767     +0.976396     +1.014447     +1.121264     
356.219   +0.982852     +0.989157     +1.006691     +1.055237     
358.018   +0.995280     +0.997020     +1.001848     +1.015101     
359.817   +0.999960     +0.999975     +1.000016     +1.000128     
SHAR_EOF
fi # end of overwriting check
if test -f 'res2'
then
	echo shar: will not over-write existing file "'res2'"
else
cat << "SHAR_EOF" > 'res2'

 typ= order n= Circumferential(c) or Radial(r):  kind: First(f) Second(d):  Standard(s) or Modified(m): 

   3 

0.125   +569.918400    +267.818838    +58.537904    +2.411248    
0.138   +427.663675    +200.753224    +43.693420    +1.735272    
0.150   +328.969182    +154.242530    +33.415947    +1.276341    
0.163   +258.367936    +120.985614    +26.080824    +0.955820    
0.175   +206.540818    +96.583969    +20.709788    +0.726635    
0.188   +167.644617    +78.280048    +16.689873    +0.559476    
0.200   +137.888921    +64.285330    +13.623755    +0.435482    
0.213   +114.741769    +53.405287    +11.246236    +0.342164    
0.225   +96.468020    +44.821452    +9.375731    +0.271045    
0.238   +81.851530    +37.960272    +7.885086    +0.216248    
0.250   +70.022740    +32.411725    +6.683470    +0.173620    
0.263   +60.348718    +27.877428    +5.704839    +0.140175    
0.275   +52.361340    +24.136738    +4.900399    +0.113738    
0.288   +45.709328    +21.024121    +4.233576    +0.092700    
0.300   +40.125602    +18.413752    +3.676596    +0.075856    
0.313   +35.404651    +16.208833    +3.208116    +0.062299    
0.325   +31.386597    +14.334081    +2.811554    +0.051333    
0.338   +27.945800    +12.730349    +2.473896    +0.042424    
0.350   +24.982607    +11.350738    +2.184838    +0.035157    
0.363   +22.417293    +10.157737    +1.936145    +0.029208    
0.375   +20.185561    +9.121105    +1.721190    +0.024321    
0.388   +18.235177    +8.216280    +1.534597    +0.020295    
0.400   +16.523403    +7.423174    +1.371975    +0.016968    
0.413   +15.015053    +6.725252    +1.229716    +0.014213    
0.425   +13.680974    +6.108821    +1.104836    +0.011925    
0.438   +12.496878    +5.562475    +0.994855    +0.010021    
0.450   +11.442420    +5.076666    +0.897699    +0.008433    
0.463   +10.500472    +4.643356    +0.811625    +0.007106    
0.475   +9.656544    +4.255749    +0.735164    +0.005996    
0.488   +8.898325    +3.908074    +0.667069    +0.005065    
0.500   +8.215308    +3.595407    +0.606280    +0.004283    
0.513   +7.598490    +3.313531    +0.551890    +0.003626    
0.525   +7.040127    +3.058820    +0.503121    +0.003073    
0.538   +6.533535    +2.828145    +0.459304    +0.002607    
0.550   +6.072922    +2.618798    +0.419860    +0.002213    
0.563   +5.653253    +2.428425    +0.384288    +0.001880    
0.575   +5.270141    +2.254976    +0.352153    +0.001599    
0.588   +4.919746    +2.096656    +0.323074    +0.001361    
0.600   +4.598701    +1.951896    +0.296721    +0.001159    
0.612   +4.304045    +1.819312    +0.272801    +0.000987    
0.625   +4.033165    +1.697687    +0.251060    +0.000842    
0.637   +3.783752    +1.585946    +0.231271    +0.000718    
0.650   +3.553760    +1.483135    +0.213238    +0.000613    
0.662   +3.341369    +1.388407    +0.196782    +0.000524    
0.675   +3.144963    +1.301011    +0.181750    +0.000448    
0.687   +2.963094    +1.220274    +0.168001    +0.000383    
0.700   +2.794473    +1.145596    +0.155413    +0.000327    
0.712   +2.637940    +1.076441    +0.143876    +0.000280    
0.725   +2.492457    +1.012326    +0.133291    +0.000240    
0.737   +2.357090    +0.952819    +0.123571    +0.000206    
0.750   +2.230996    +0.897529    +0.114636    +0.000176    
0.762   +2.113414    +0.846104    +0.106416    +0.000151    
0.775   +2.003658    +0.798228    +0.098848    +0.000130    
0.787   +1.901105    +0.753612    +0.091873    +0.000111    
0.800   +1.805190    +0.711996    +0.085441    +0.000095    
0.812   +1.715401    +0.673144    +0.079504    +0.000082    
0.825   +1.631270    +0.636840    +0.074020    +0.000070    
0.837   +1.552372    +0.602889    +0.068952    +0.000060    
0.850   +1.478319    +0.571112    +0.064264    +0.000052    
0.862   +1.408758    +0.541347    +0.059924    +0.000045    
0.875   +1.343362    +0.513445    +0.055906    +0.000038    
0.887   +1.281836    +0.487270    +0.052181    +0.000033    
0.900   +1.223906    +0.462697    +0.048728    +0.000028    
0.912   +1.169323    +0.439612    +0.045523    +0.000024    
0.925   +1.117856    +0.417910    +0.042548    +0.000021    
0.937   +1.069293    +0.397494    +0.039785    +0.000018    
0.950   +1.023440    +0.378275    +0.037216    +0.000016    
0.962   +0.980116    +0.360173    +0.034828    +0.000013    
0.975   +0.939155    +0.343110    +0.032606    +0.000012    
0.987   +0.900404    +0.327019    +0.030537    +0.000010    
1.000   +0.863720    +0.311833    +0.028611    +0.000009    
1.012   +0.828971    +0.297495    +0.026816    +0.000007    
1.025   +0.796036    +0.283948    +0.025143    +0.000006    
1.037   +0.764802    +0.271142    +0.023583    +0.000006    
1.050   +0.735162    +0.259029    +0.022127    +0.000005    
1.062   +0.707019    +0.247567    +0.020768    +0.000004    
1.075   +0.680283    +0.236713    +0.019499    +0.000004    
1.087   +0.654868    +0.226430    +0.018313    +0.000003    
1.100   +0.630696    +0.216683    +0.017205    +0.000003    
1.112   +0.607694    +0.207440    +0.016169    +0.000002    
1.125   +0.585792    +0.198669    +0.015200    +0.000002    
1.137   +0.564926    +0.190343    +0.014293    +0.000002    
1.150   +0.545037    +0.182434    +0.013444    +0.000001    
1.162   +0.526068    +0.174918    +0.012649    +0.000001    
1.175   +0.507968    +0.167773    +0.011905    +0.000001    
1.187   +0.490687    +0.160976    +0.011207    +0.000001    
1.200   +0.474180    +0.154508    +0.010552    +0.000001    
1.212   +0.458403    +0.148349    +0.009939    +0.000001    
1.225   +0.443316    +0.142483    +0.009363    +0.000001    
1.237   +0.428882    +0.136892    +0.008823    +0.000001    
1.250   +0.415064    +0.131561    +0.008316    +0.000000    
1.262   +0.401831    +0.126476    +0.007840    +0.000000    
1.275   +0.389149    +0.121623    +0.007393    +0.000000    
1.287   +0.376991    +0.116990    +0.006973    +0.000000    
1.300   +0.365329    +0.112564    +0.006578    +0.000000    
1.312   +0.354136    +0.108335    +0.006207    +0.000000    
1.325   +0.343388    +0.104292    +0.005857    +0.000000    
1.337   +0.333063    +0.100425    +0.005529    +0.000000    
1.350   +0.323138    +0.096725    +0.005220    +0.000000    
1.362   +0.313594    +0.093183    +0.004929    +0.000000    
1.375   +0.304411    +0.089792    +0.004655    +0.000000    
1.387   +0.295571    +0.086543    +0.004397    +0.000000    
1.400   +0.287058    +0.083429    +0.004154    +0.000000    
1.412   +0.278855    +0.080444    +0.003925    +0.000000    
1.425   +0.270947    +0.077581    +0.003709    +0.000000    
1.437   +0.263320    +0.074833    +0.003505    +0.000000    
1.450   +0.255961    +0.072197    +0.003314    +0.000000    
1.462   +0.248857    +0.069665    +0.003133    +0.000000    
1.475   +0.241996    +0.067233    +0.002962    +0.000000    
1.487   +0.235366    +0.064896    +0.002801    +0.000000    
1.500   +0.228958    +0.062650    +0.002649    +0.000000    
1.512   +0.222760    +0.060490    +0.002505    +0.000000    
1.525   +0.216765    +0.058413    +0.002370    +0.000000    
1.537   +0.210961    +0.056415    +0.002242    +0.000000    
1.550   +0.205342    +0.054492    +0.002121    +0.000000    
1.562   +0.199899    +0.052640    +0.002007    +0.000000    
1.575   +0.194624    +0.050857    +0.001899    +0.000000    
1.587   +0.189510    +0.049139    +0.001798    +0.000000    
1.600   +0.184550    +0.047484    +0.001701    +0.000000    
1.612   +0.179739    +0.045889    +0.001610    +0.000000    
1.625   +0.175069    +0.044351    +0.001524    +0.000000    
1.637   +0.170535    +0.042868    +0.001443    +0.000000    
1.650   +0.166132    +0.041438    +0.001366    +0.000000    
1.662   +0.161854    +0.040058    +0.001293    +0.000000    
1.675   +0.157697    +0.038726    +0.001224    +0.000000    
1.687   +0.153655    +0.037441    +0.001159    +0.000000    
1.700   +0.149724    +0.036200    +0.001098    +0.000000    
1.712   +0.145900    +0.035002    +0.001040    +0.000000    
1.725   +0.142179    +0.033844    +0.000984    +0.000000    
1.737   +0.138558    +0.032726    +0.000932    +0.000000    
1.750   +0.135032    +0.031646    +0.000883    +0.000000    
1.762   +0.131598    +0.030602    +0.000836    +0.000000    
1.775   +0.128252    +0.029592    +0.000792    +0.000000    
1.787   +0.124993    +0.028617    +0.000750    +0.000000    
1.800   +0.121816    +0.027674    +0.000710    +0.000000    
1.812   +0.118720    +0.026761    +0.000673    +0.000000    
1.825   +0.115700    +0.025879    +0.000637    +0.000000    
1.837   +0.112756    +0.025026    +0.000603    +0.000000    
1.850   +0.109883    +0.024200    +0.000571    +0.000000    
1.862   +0.107081    +0.023401    +0.000541    +0.000000    
1.875   +0.104346    +0.022628    +0.000513    +0.000000    
1.887   +0.101678    +0.021880    +0.000485    +0.000000    
1.900   +0.099072    +0.021156    +0.000460    +0.000000    
1.912   +0.096529    +0.020456    +0.000435    +0.000000    
1.925   +0.094046    +0.019777    +0.000412    +0.000000    
1.937   +0.091621    +0.019120    +0.000391    +0.000000    
1.950   +0.089253    +0.018485    +0.000370    +0.000000    
1.962   +0.086940    +0.017869    +0.000350    +0.000000    
1.975   +0.084680    +0.017273    +0.000332    +0.000000    
1.987   +0.082473    +0.016696    +0.000314    +0.000000    
2.000   +0.080317    +0.016137    +0.000297    +0.000000    
2.012   +0.078210    +0.015595    +0.000282    +0.000000    
2.025   +0.076151    +0.015071    +0.000267    +0.000000    
2.037   +0.074140    +0.014563    +0.000252    +0.000000    
2.050   +0.072174    +0.014072    +0.000239    +0.000000    
2.062   +0.070254    +0.013596    +0.000226    +0.000000    
2.075   +0.068377    +0.013135    +0.000214    +0.000000    
2.087   +0.066544    +0.012688    +0.000203    +0.000000    
2.100   +0.064752    +0.012256    +0.000192    +0.000000    
2.112   +0.063002    +0.011837    +0.000182    +0.000000    
2.125   +0.061291    +0.011432    +0.000172    +0.000000    
2.137   +0.059620    +0.011039    +0.000163    +0.000000    
2.150   +0.057988    +0.010659    +0.000154    +0.000000    
2.162   +0.056393    +0.010291    +0.000146    +0.000000    
2.175   +0.054836    +0.009935    +0.000138    +0.000000    
2.187   +0.053315    +0.009590    +0.000131    +0.000000    
2.200   +0.051829    +0.009256    +0.000123    +0.000000    
2.212   +0.050378    +0.008933    +0.000117    +0.000000    
2.225   +0.048962    +0.008620    +0.000111    +0.000000    
2.237   +0.047579    +0.008318    +0.000105    +0.000000    
2.250   +0.046230    +0.008025    +0.000099    +0.000000    
2.262   +0.044912    +0.007742    +0.000094    +0.000000    
2.275   +0.043627    +0.007468    +0.000089    +0.000000    
2.287   +0.042373    +0.007203    +0.000084    +0.000000    
2.300   +0.041150    +0.006946    +0.000079    +0.000000    
2.312   +0.039956    +0.006698    +0.000075    +0.000000    
2.325   +0.038793    +0.006459    +0.000071    +0.000000    
2.337   +0.037658    +0.006227    +0.000067    +0.000000    
2.350   +0.036552    +0.006003    +0.000063    +0.000000    
2.363   +0.035474    +0.005786    +0.000060    +0.000000    
2.375   +0.034423    +0.005577    +0.000057    +0.000000    
2.388   +0.033400    +0.005375    +0.000053    +0.000000    
2.400   +0.032403    +0.005179    +0.000051    +0.000000    
2.413   +0.031431    +0.004991    +0.000048    +0.000000    
2.425   +0.030486    +0.004808    +0.000045    +0.000000    
2.438   +0.029565    +0.004632    +0.000043    +0.000000    
2.450   +0.028668    +0.004462    +0.000040    +0.000000    
2.463   +0.027796    +0.004298    +0.000038    +0.000000    
2.475   +0.026947    +0.004139    +0.000036    +0.000000    
2.488   +0.026121    +0.003986    +0.000034    +0.000000    
2.500   +0.025318    +0.003839    +0.000032    +0.000000    
2.513   +0.024537    +0.003696    +0.000030    +0.000000    
2.525   +0.023777    +0.003558    +0.000029    +0.000000    
2.538   +0.023039    +0.003426    +0.000027    +0.000000    
2.550   +0.022321    +0.003298    +0.000026    +0.000000    
2.563   +0.021623    +0.003174    +0.000024    +0.000000    
2.575   +0.020945    +0.003055    +0.000023    +0.000000    
2.588   +0.020287    +0.002940    +0.000022    +0.000000    
2.600   +0.019647    +0.002830    +0.000020    +0.000000    
2.613   +0.019026    +0.002723    +0.000019    +0.000000    
2.625   +0.018422    +0.002620    +0.000018    +0.000000    
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'bsslr.cpp'
then
	echo shar: will not over-write existing file "'bsslr.cpp'"
else
cat << "SHAR_EOF" > 'bsslr.cpp'
#include   <iostream>
#include   <stdlib.h>
#include   <cmath>

/*
    Standard Bessel functions of first and second kind.
    Computation is based on routines from Numerical Recipes in C with
    some modifications.
    In most cases accuracy is more than 14 decimal places.
    The routines are specifically designed to compute arrays of
    Jn,J'n,Yn and Y'n

*/

inline int IMAX(int a,int b)
	{if (a>b) return a; else return b;}
inline double SIGN(double a, double b)
	      { if (b>=0.0) return fabs(a); else return -fabs(a);}
inline double abs(double z) {return fabs(z);}
const double pi=3.141592653589793;

       /* Routines in this file*/
       double lnGamma(double xx);
       double Gamma(double n);
       double Gamma(int  n);

       double BesselJn(double order, double z);
       double BesselJDn(double order, double z);
       double BesselJ(double z, int M, double bJn[]);
       double BesselJ(double z, int M, double bJn[], double JDn[]);

	double BesselYn(double order, double z);
	double BesselYDn(double order, double z);
	double BesselY(double z, int M, double bYn[]);
	double BesselY(double z, int M, double bYn[], double bYDn[]);

	void bessjy(double x, double xnu, double *Jv, double *Yv, double *JDv, double *YDv); /* fractional order as well as integer */

	double bessj(int n, double x);
	double bessjD(int n, double x);
	double bessy(int n, double x);
	double bessyD(int n, double x);

	/* Arrays of fractional order */
	double BesselJ(double z,double v, int M, double pJn[],double mJn[]);

/*------------- Gamma Function  ---------------------*/
double lnGamma(double xx)
{
	double x,y,tmp,ser;
	static double cof[6]={76.18009172947146,-86.50532032941677,
		24.01409824083091,-1.231739572450155,
		0.1208650973866179e-2,-0.5395239384953e-5};
	int j;

	y=x=xx;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+log(2.5066282746310005*ser/x);
}

double Gamma(double n)  /*   G(n)=(n-1)!  */
{
  int i,nend=20;
  double z,apx,gn;

  if (n==1) return 1;
  if (n+int(abs(n))==0) return 1e300;
  nend=nend+int(abs(n));
  z=n+nend;
  apx=exp((z-.5)*log(z)-z+.5*log(2*pi)+1/(12*z)-pow(z,-3)/360+pow(z,-5)/1260-pow(z,-7)/1680);
  gn=apx;
  for (i=nend;i>=1; i--)
  { z=z-1;
    gn=gn/z;
   }
   return gn;
}

double  Gamma(int n)  /*   G(n)=(n-1)!  */
{
  double gn,i;

  if (n<=0) return 1e300;
  gn=1;
  for (i=1;i<n; i++)   gn=gn*i;

   return gn;
}
/*===============================================================*/
/*--------------- Bessel functions Jn Real----------------------*/
double BesselJn(double order , double z)
{
   double Jn,Yn,JDn,YDn,acc;

   if (z==0)   { if (order==0) return 1; else return 0;}

	bessjy(z,order,&Jn,&Yn,&JDn,&YDn);
	acc=fabs((Jn*YDn-JDn*Yn)*0.5*z*pi-1.0);
    if(acc>1.0e-9) {std::cout <<"\n BesselJn: Compuation is not accurate: acc="<<acc;}

  return  Jn;

}
/*------------------------------------------------*/
double BesselJDn(double order , double z)
{
   double Jn,Yn,JDn,YDn,acc;

   if (z==0)
   { if (order==1) return 0.5; else return 0;}

	bessjy(z,order,&Jn,&Yn,&JDn,&YDn);
	acc=fabs((Jn*YDn-JDn*Yn)*0.5*z*pi-1.0);
    if(acc>1.0e-9) {std::cout <<"\n BesselJDn: Compuation is not accurate: acc="<<acc;}

   return  JDn;

}
/*------------------ Arrays of Bessel Real -----------------*/
double BesselJ(double z, int M, double bJn[])
{       /* Array size is M, i.e bJn[0],...,bJn[M-1]  */
  int i;

  if (z==0)
   { bJn[0]=1;
     for (i=1;i<M;i++)
      bJn[i]=0;
      return 1e-20;
   }
	int j,SJ;
	double az,bj,bjm,bjp,toz,fact;
	double ACR=72.0,sign,Sum;
	double J0,Y0,JD0,YD0,acc=0,J1;

	az=fabs(z);
	toz=2.0/az;

	if (az > (double) M) {
		bessjy(az,0,&J0,&Y0,&JD0,&YD0);
		acc=fabs((J0*YD0-JD0*Y0)*0.5*az*pi-1.0);
		J1=-JD0;
		bJn[0]=J0;
		bJn[1]=J1;
		for (j=1;j<M-1;j++) {
			bJn[j+1]=j*toz*bJn[j]-bJn[j-1];
		}
	} else {
		for(i=0;i<M;i++)   bJn[i]=0.0;
		Sum=0.0;
		bjp=0.0;
		bj=1.0;
		SJ=int(2*((M+int(sqrt(ACR*M))/2.0)));
		for (j=SJ;j>0;j--) {
			bjm=j*toz*bj-bjp;
			bjp=bj;
			bj=bjm;
			if (fabs(bj) > 1.0e10) {
			     for(i=0;i<M;i++)   bJn[i] *=1.0e-10;
				Sum *=1.0e-10;
				bj *= 1.0e-10;
				bjp *= 1.0e-10;
			}
			if(j%2) Sum +=bj;
			if (j<M) bJn[j]=bjp;
		}
	     if(az<1)
		   {
		     bJn[0]=bj;
		     Sum=2.0*Sum-bJn[0];
		     J0=bJn[0]/Sum;
		   }
		else
		  {
		    bessjy(az,0,&J0,&Y0,&JD0,&YD0);
		    acc=fabs((J0*YD0-JD0*Y0)*0.5*az*pi-1.0);
		  }

		fact=J0/bj;
		bJn[0]=J0;
		sign=1.0;
		for(j=1;j<M;j++)
		{  if(z<0.0) sign=-sign;
		   bJn[j]=sign*bJn[j]*fact;
		}
	}
   if(acc>1.0e-9) {std::cout <<"\n BesselJ: Compuation is not accurate: acc="<<acc;}
    return acc;
}
/*-------------------------------------------------------*/
double BesselJ(double z, int M, double bJn[], double JDn[])
{
   int n;
   double acc;

   acc=BesselJ(z,M,bJn);
   JDn[0]=-bJn[1];
   for(n=1;n<M-1;n++)
    { JDn[n]=0.5*(bJn[n-1]-bJn[n+1]);}
   if(z==0) JDn[M-1]=0; else JDn[M-1]=bJn[M-2]-(M-1)*bJn[M-1]/z;

	return acc;
}

/*====================  Bessel Y   =================================*/
/*------------- Bessel Y real -------------------------------------*/
double BesselYn(double order, double z)
{
   double Jn,Yn,JDn,YDn,acc;
  if (abs(z)==0) return 1e200;

	bessjy(z,order,&Jn,&Yn,&JDn,&YDn);
	acc=fabs((Jn*YDn-JDn*Yn)*0.5*z*pi-1.0);
   if(acc>1.0e-9) {std::cout <<"\n BesselYn: Compuation is not accurate: acc="<<acc;}
    return Yn;

}
double BesselYDn(double order, double z)
{
   double Jn,Yn,JDn,YDn,acc;
  if (abs(z)==0) return 1e200;

	bessjy(z,order,&Jn,&Yn,&JDn,&YDn);
	acc=fabs((Jn*YDn-JDn*Yn)*0.5*z*pi-1.0);
    if(acc>1.0e-9) {std::cout <<"\n BesselYDn: Compuation is not accurate: acc="<<acc;}

    return YDn;

}
/*------------------ Bessel Arrays Y real -------------------------*/
double BesselY(double z, int M, double bYn[])
{

 if (z==0)
     {
      std::cout <<"BesselY d: zero argument i.e infinit result can not containue";
      return 1e9;
      }

	double by,bym,byp,toz,acc;
	double Jn,Yn,JDn,YDn;
	int j,i;

	toz=2.0/z;
	bessjy(z,0,&Jn,&Yn,&JDn,&YDn);
	acc=fabs((Jn*YDn-JDn*Yn)*0.5*z*pi-1.0);
	by=-YDn; // for n=1
	bym=Yn; // for n=0
	bYn[0]=bym;
	bYn[1]=by;
	for (j=1;j<M-1;j++) {
		byp=j*toz*by-bym;
		bym=by;
		by=byp;
		bYn[j+1]=by;
	  if(j>4 && fabs(by)>1e260) {break;}
//	       {std::cout<<"BesselY: Yn array is larger than largest machine number z="<<z;exit(0);}
	}
	if(j<M-1)
	      for(i=j+2;i<M;i++) bYn[i]=1e260;

	if(acc>1.0e-9) {std::cout <<"\n BesselY: Compuation is not accurate: acc="<<acc;}
	return acc;

}
double BesselY(double z, int M, double bYn[], double bYDn[])
{
   int n;
   double acc;

    if (z==0)
     {
      std::cout <<"BesselY d: zero argument i.e infinit result can not containue";
      return 1e9;
      }

   acc=BesselY(z,M,bYn);
   bYDn[0]=-bYn[1];
   for(n=1;n<M-1;n++)
    { bYDn[n]=0.5*(bYn[n-1]-bYn[n+1]);}
  bYDn[M-1]=bYn[M-2]-(M-1.0)*bYn[M-1]/z;

	return acc;
}
/*================================================================*/
/*-----------------------------------------------------------*/

double bessj(int n, double z)
{
	int j,jsum,m;
	double az,bj,bjm,bjp,sum,toz,ans,sign=1.0;
	double ACR=400.0;
	double J0,Y0,JD0,YD0,acc,J1;

	if (z == 0.0)	if(n==0) return 1; else return 0.0;
	az=fabs(z);
	if(z<0) sign=pow(-1,n);
	else if (az > (double) n) {
		bessjy(az,0,&J0,&Y0,&JD0,&YD0);
		acc=fabs((J0*YD0-JD0*Y0)*0.5*z*pi-1.0);
		J1=-JD0;
		if(n==0) return J0;
		if(n==1) return sign*J1;
		toz=2.0/az;
		bjm=J0;
		bj=J1;
		for (j=1;j<n;j++) {
			bjp=j*toz*bj-bjm;
			bjm=bj;
			bj=bjp;
		}
		ans=sign*bj;
	} else {
		toz=2.0/az;
		m=int(2*((n+(int) sqrt(ACR*n))/2.0));
		jsum=0;
		bjp=ans=sum=0.0;
		bj=1.0;
		for (j=m;j>0;j--) {
			bjm=j*toz*bj-bjp;
			bjp=bj;
			bj=bjm;
			if (fabs(bj) > 1.0e10) {
				bj *= 1.0e-10;
				bjp *= 1.0e-10;
				ans *= 1.0e-10;
				sum *= 1.0e-10;
			}
			if (jsum) sum += bj;
			jsum=!jsum;
			if (j == n) ans=bjp;
		}
		sum=2.0*sum-bj;
		ans /= sum;
		ans=ans*sign;
	}
	if(acc>1.0e-9) {std::cout <<"\n bessj: Compuation is not accurate: acc="<<acc;}

	return ans;
}

/*====================================================*/
double bessy(int n, double z)
{
	int j;
	double by,bym,byp,toz;
	double J0,Y0,JD0,YD0,acc,Y1;

	bessjy(z,0,&J0,&Y0,&JD0,&YD0);
	acc=fabs((J0*YD0-JD0*Y0)*0.5*z*pi-1.0);
	Y1=-YD0;
	if(n==0) return Y0;
	if(n==1) return Y1;

	toz=2.0/z;
	by=Y1;
	bym=Y0;
	for (j=1;j<n;j++) {
		byp=j*toz*by-bym;
		bym=by;
		by=byp;
	}
	if(acc>1.0e-9) {std::cout <<"\n bessy: Compuation is not accurate: acc="<<acc;}
	return by;
}
double bessyD(int n, double z)
{
	int j;
	double by,bym,byp,toz;
	double J0,Y0,JD0,YD0,acc,Y1;

	bessjy(z,0,&J0,&Y0,&JD0,&YD0);
	acc=fabs((J0*YD0-JD0*Y0)*0.5*z*pi-1.0);
	Y1=-YD0;
	if(n==0) return -YD0;

	toz=2.0/z;
	by=Y1;
	bym=Y0;
	for (j=1;j<=n;j++) {
		byp=j*toz*by-bym;
		bym=by;
		by=byp;
	}
	if(acc>1.0e-9) {std::cout <<"\n bessyD: Compuation is not accurate: acc="<<acc;}

	return bym*n/z-byp;
}
/*-------------------------------------------------------------------*/
/*         Fractional Order Bessel Function                          */
/*-------------------------------------------------------------------*/
void bessjy(double x, double xnu, double *Jv, double *Yv, double *JDv, double *YDv)
{
	void beschB(double x, double *gam1, double *gam2, double *gampl,
		double *gammi);

	int i,isign,l,nl;
	double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
		fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,Jvl,
		Jvl1,Jvmu,JDv1,JDvl,Jvtemp,Yv1,Yvmu,Yvmup,Yvtemp,sum,sum1,
		temp,w,x2,xi,xi2,xmu,xmu2;
	double EPS=1.0e-16,FPMIN=1.0e-30;
	int MAXIT=10000,XMIN=2;

	if (x <= 0.0 || xnu < 0.0) {std::cout<<"bad arguments in bessjy";exit(0);}
	nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5)));
	xmu=xnu-nl;
	xmu2=xmu*xmu;
	xi=1.0/x;
	xi2=2.0*xi;
	w=xi2/pi;
	isign=1;
	h=xnu*xi;
	if (h < FPMIN) h=FPMIN;
	b=xi2*xnu;
	d=0.0;
	c=h;
	for (i=1;i<=MAXIT;i++) {
		b += xi2;
		d=b-d;
		if (fabs(d) < FPMIN) d=FPMIN;
		c=b-1.0/c;
		if (fabs(c) < FPMIN) c=FPMIN;
		d=1.0/d;
		del=c*d;
		h=del*h;
		if (d < 0.0) isign = -isign;
		if (fabs(del-1.0) < EPS) break;
	}
	if (i > MAXIT) {std::cout<<"x too large in bessjy; tYv asymptotic expansion";exit(0);}
	Jvl=isign*FPMIN;
	JDvl=h*Jvl;
	Jvl1=Jvl;
	JDv1=JDvl;
	fact=xnu*xi;
	for (l=nl;l>=1;l--) {
		Jvtemp=fact*Jvl+JDvl;
		fact -= xi;
		JDvl=fact*Jvtemp-Jvl;
		Jvl=Jvtemp;
	}
	if (Jvl == 0.0) Jvl=EPS;
	f=JDvl/Jvl;
	if (x < XMIN) {
		x2=0.5*x;
		pimu=pi*xmu;
		fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
		d = -log(x2);
		e=xmu*d;
		fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
		beschB(xmu,&gam1,&gam2,&gampl,&gammi);
		ff=2.0/pi*fact*(gam1*cosh(e)+gam2*fact2*d);
		e=exp(e);
		p=e/(gampl*pi);
		q=1.0/(e*pi*gammi);
		pimu2=0.5*pimu;
		fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2);
		r=pi*pimu2*fact3*fact3;
		c=1.0;
		d = -x2*x2;
		sum=ff+r*q;
		sum1=p;
		for (i=1;i<=MAXIT;i++) {
			ff=(i*ff+p+q)/(i*i-xmu2);
			c *= (d/i);
			p /= (i-xmu);
			q /= (i+xmu);
			del=c*(ff+r*q);
			sum += del;
			del1=c*p-i*del;
			sum1 += del1;
			if (fabs(del) < (1.0+fabs(sum))*EPS) break;
		}
		if (i > MAXIT) {std::cout<<"bessy series failed to converge";exit(0);}
		Yvmu = -sum;
		Yv1 = -sum1*xi2;
		Yvmup=xmu*xi*Yvmu-Yv1;
		Jvmu=w/(Yvmup-f*Yvmu);
	} else {
		a=0.25-xmu2;
		p = -0.5*xi;
		q=1.0;
		br=2.0*x;
		bi=2.0;
		fact=a*xi/(p*p+q*q);
		cr=br+q*fact;
		ci=bi+p*fact;
		den=br*br+bi*bi;
		dr=br/den;
		di = -bi/den;
		dlr=cr*dr-ci*di;
		dli=cr*di+ci*dr;
		temp=p*dlr-q*dli;
		q=p*dli+q*dlr;
		p=temp;
		for (i=2;i<=MAXIT;i++) {
			a += 2*(i-1);
			bi += 2.0;
			dr=a*dr+br;
			di=a*di+bi;
			if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN;
			fact=a/(cr*cr+ci*ci);
			cr=br+cr*fact;
			ci=bi-ci*fact;
			if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN;
			den=dr*dr+di*di;
			dr /= den;
			di /= -den;
			dlr=cr*dr-ci*di;
			dli=cr*di+ci*dr;
			temp=p*dlr-q*dli;
			q=p*dli+q*dlr;
			p=temp;
			if (fabs(dlr-1.0)+fabs(dli) < EPS) break;
		}
		if (i > MAXIT)
		     { std::cout <<"cf2 failed in bessjy"; exit(0); }
		gam=(p-f)/q;
		Jvmu=sqrt(w/((p-f)*gam+q));
		Jvmu=SIGN(Jvmu,Jvl);
		Yvmu=Jvmu*gam;
		Yvmup=Yvmu*(p+q/gam);
		Yv1=xmu*xi*Yvmu-Yvmup;
	}
	fact=Jvmu/Jvl;
	*Jv=Jvl1*fact;
	*JDv=JDv1*fact;
	for (i=1;i<=nl;i++) {
		Yvtemp=(xmu+i)*xi2*Yv1-Yvmu;
		Yvmu=Yv1;
		Yv1=Yvtemp;
	}
	*Yv=Yvmu;
	*YDv=xnu*xi*Yvmu-Yv1;
}
void beschB(double x, double *gam1, double *gam2, double *gampl, double *gammi)
{
	double Chebev(double a, double b, double c[], int m, double x);
	double xx;
	int NUSE1=7,NUSE2=8;
	static double c1[] = {
		-1.142022680371172e0,6.516511267076e-3,
		3.08709017308e-4,-3.470626964e-6,6.943764e-9,
		3.6780e-11,-1.36e-13};
	static double c2[] = {
		1.843740587300906e0,-0.076852840844786e0,
		1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
		2.42310e-10,-1.70e-13,-1.0e-15};

	xx=8.0*x*x-1.0;
	*gam1=Chebev(-1.0,1.0,c1,NUSE1,xx);
	*gam2=Chebev(-1.0,1.0,c2,NUSE2,xx);
	*gampl= *gam2-x*(*gam1);
	*gammi= *gam2+x*(*gam1);
}
double Chebev(double a, double b, double c[], int m, double x)
{
	double d=0.0,dd=0.0,sv,y,y2;
	int j;

	if ((x-a)*(x-b) > 0.0) {std::cout<<"x not in range in routine Chebev";exit(0);}
	y2=2.0*(y=(2.0*x-a-b)/(b-a));
	for (j=m-1;j>=1;j--) {
		sv=d;
		d=y2*d-dd+c[j];
		dd=sv;
	}
	return y*d-dd+0.5*c[0];
}
/*-------------------------------------------------------------------------*/
/*  Takes care of negative orders  */
void nbessjy(double x, double xnu, double *Jv, double *Yv, double *JDv, double *YDv)
{

	double	Jvl1,Jvtemp,Yv1,Yvtemp,xmu;
	int sign=1;

	xmu=xnu;
	if (xmu<0)  {xmu=-xmu;sign=-1;}
	bessjy(x,xmu,Jv,Yv,JDv,YDv);

	if(sign<0)  /* use reflection formula */
	{
	 Jvl1= *Jv*cos(xmu*pi) - *Yv*sin(xmu*pi);
	 Yv1= *Jv*sin(xmu*pi) + *Yv*cos(xmu*pi);
	 Jvtemp= *JDv*cos(xmu*pi) - *YDv*sin(xmu*pi);
	 Yvtemp= *JDv*sin(xmu*pi) + *YDv*cos(xmu*pi);
	 *Jv=Jvl1;
	 *Yv=Yv1;
	 *JDv=Jvtemp;
	 *YDv=Yvtemp;
	}
}
/*===================================================================*/
SHAR_EOF
fi # end of overwriting check
if test -f 'bsslr.h'
then
	echo shar: will not over-write existing file "'bsslr.h'"
else
cat << "SHAR_EOF" > 'bsslr.h'
       double Gamma(double n);
       double Gamma(int  n);

       double BesselJn(double order, double z);
       double BesselJDn(double order, double z);
       double BesselJ(double z, int M, double bJn[]);
       double BesselJ(double z, int M, double bJn[], double JDn[]);

	double BesselYn(double order, double z);
	double BesselYDn(double order, double z);
	double BesselY(double z, int M, double Yn[]);
	double BesselY(double z, int M, double Yn[], double YDn[]);
SHAR_EOF
fi # end of overwriting check
if test -f 'mathur.cpp'
then
	echo shar: will not over-write existing file "'mathur.cpp'"
else
cat << "SHAR_EOF" > 'mathur.cpp'
/* standard include files*/
#include <iostream>
#include <cstdlib>
#include <cmath>

#define CDIMS 60  /* Important: CIDMs must be greater than the order n, i.e. CDIMs=n+15 */
/*
  The Author would be happy to help in any way he can
  in the use of these routins.

  Contact Address:
  CERI, KACST, P.O. Box 6086,
  Riyadh 11442,
  Saudi Arabia,
  Fax:966+1+4813764
  Email: alhargan@kacst.edu.sa

  Legal Matters:
  The program is here as is, no guarantees are given, stated or implied,
  you may use it at your own risk.

*/
/*
  For theoretical derivation see:
  The companion paper
  or
  F.A. Alhargan, "A complete method for the computations of Mathieu
  characteristic numbers of integer orders," SIAM Review, vol-38, No. 2,
  pp. 239-255, June 1996.

*/
/*
  Mathieu functions of real argument
  last modified 20-12-1998

*/
/*
  To be able to run any of these routines you need the following programs:
  bsslr.cpp
  mcnr.cpp

*/

#define CDIMs 60    /* this should be higher than the order n, i.e CDIMs=n+15 */


const  double pi=3.141592653589793;

inline double abs(double z)
     {return  fabs(z);}
inline double acosh(double z)
       {return log(z+sqrt(z*z-1));}
inline double asinh(double z)
       {return log(z+sqrt(z*z+1));}
inline double atanh(double z)
       {return 0.5*log((1+z)/(1-z));}
inline int signf(int m)     /* sign function */
{
  if(m%2==0) return 1;
   else return -1;
}


#define h_type double
	/*
	  For complex h replace double by complex, routins in this file
	  do not need modifing. However, you will not be able to run the
	  complex  part unless you have the complex routins for computing
	  Mathieu  Charactristic Numbers (MCNs) mcnc.cpp.
	  The one I have just developed  has limited range for h complex,
	  I am working in its extension.
	  Though it will not be  available  for sometime.
	*/

     /* routines required from file mcnr.cpp  */
   h_type Estimatmcn(char typ,int n, h_type h);
   h_type MCNRoot(h_type a, char typ, int n, h_type h, double accr, double &acc);
   double Coefficients(char typ,int n, h_type h, h_type Bm[], int CDIM, h_type &mcnr, h_type &Norm);
   int CoefficientSec(char typ, int n, h_type h, h_type mcn, h_type Bm[], int CDIM, h_type Dm[], h_type &Norm);

   /* routines required from file bsslr.cpp */
       double BesselJ(double z, int M, double bJn[]);
       double BesselJ(double z, int M, double bJn[], double JDn[]);

       double BesselY(double z, int M, double Yn[]);
       double BesselY(double z, int M, double Yn[], double YDn[]);

   /*   Routines available in this file  */
   h_type MathuZn(char typ,int n, h_type h, double u, short kind);
   h_type MathuZDn(char typ,int n, h_type h, double u, short kind);
   h_type MathuSn(char typ,int n, h_type h, double v, short kind);
   h_type MathuSDn(char typ,int n, h_type h, double v, short kind);

   h_type  MathuZn(char typ, int n, h_type Bm[], h_type Jms[], h_type Zm[], int CDIM);
   h_type  MathuZDn(char typ, int n, h_type h, double u, h_type Bm[], h_type Jms[], h_type JDms[], h_type Zm[], h_type ZDm[], int CDIM);

   h_type  MathuSn(char typ, int n, double v, h_type Bm[], int CDIM);
   h_type  MathuSDn(char typ, int n, double v, h_type Bm[], int CDIM);
   h_type  ISn(char typ,int ni, double vi, double ti, double ui, h_type Bm[], int CDIM);
   h_type  MathuFn(char typ, int n, double v, h_type Dm[], h_type Norm, h_type Bm[], int CDIM);
   h_type  MathuFDn(char typ, int n, double v, h_type Dm[], h_type Norm, h_type Bm[], int CDIM);

     /* The following do not converge well for all u, better to use MathuZn */
   h_type  MathuJn(char typ, int n, double u, h_type Bm[], h_type Jm[], int CDIM);
   h_type  MathuJDn(char typ, int n, h_type h, double u, h_type Bm[], h_type Jm[], h_type JDm[], int CDIM);
   h_type  MathuYn(char typ, int n, h_type Bm[], h_type Jms[], h_type Ym[], int CDIM);
   h_type  MathuYDn(char typ, int n, h_type h, double u, h_type Bm[], h_type Jms[], h_type JDms[], h_type Ym[], h_type YDm[], int CDIM);


/* ====================================================================*/
/* ================== Mathieu Functions ===============================*/
/*
  Important Convensions
1. CDIM is the number of terms summed in the series, which also determins
   the size of the arrays.
2. The value p is used to identify if n is even (p=0) or odd (p=1).
3. D in any function means differential e.g JDn(x)= d(Jn(x))/dx=J'n(x)
4. The series used are that which are convergent in all the z-plane.
   This is MahtuZm.
5. The Coefficients of the Even functions is assumed to be stored in  one array
   Be() for both even and odd orders, similarly for odd functions in Bo().
6. The Coefficients Be() and Bo() have to be calculated before calling the
   routines  for Mathieu functions.
7. In the case of Jen and Jon the Bessel functions have to be calculated first,
   also for Yen and Yon.
8. Warning are not serious, but care needs to be taken as to the validity of
   the results, in most cases the results are correct.
9. In the case of a serious error the programs haults with a message

10.  Bm stands for Bem or Bom     1st kind coefficients
     Dm stands for Dem or Dom     2nd kind coefficients

11. All arrays are of size 2*CDIM, i.e Bm[0]...Bm[2*CDIM-1]
    some arrays are of size CDIM, e.g. Zm[0]...Zm[CDIM-1]

*/
/*------------------------------------------------------------------*/
/*---------------- Easy Functions -----------------------------------*/
/*
  When computing Mathieu functions, one usually would be requiring several
  different functions of the same order, in this case it is best to use
  the functions which take the coefficients and arrays of Bessel functions
  as inputs.
  On the other hand, one may want to compute a single function of one order
  in this case, 'The Easy Functions' coming shortly may be better and easily
  used.
*/
h_type MathuZn(char typ,int n, h_type h, double u, short kind)
{
/*
   inputs:
   typ ='e' or 'o'
   n  = order
   h  = parameter
   u  = argument 
   kind = 1 or 2

   output:
   The function returns the value of the required  radial function


*/

  int CDIM=CDIMs-1;
  h_type Bm[2*CDIMs],Jms[2*CDIMs],Zm[2*CDIMs];
  h_type mcnr,Norm;
  h_type xp,xs;

  if(kind<1 || kind>2)
     {std::cout <<"MathuZn: The kind has not been specified correctly:STOP"; exit(0);}

  xp=0.5*h*exp(u);
  xs=0.5*h*exp(-u);
  Coefficients(typ,n,h,Bm,CDIM,mcnr,Norm);

  BesselJ(xs,2*CDIM,Jms);
  if(kind==1) BesselJ(xp,2*CDIM,Zm);
  if(kind==2) BesselY(xp,2*CDIM,Zm);

   return MathuZn(typ,n,Bm,Jms,Zm,CDIM);

}
h_type MathuZDn(char typ,int n, h_type h, double u, short kind)
{
  int CDIM=CDIMs-1;
  h_type Bm[2*CDIMs];
  h_type Jms[2*CDIMs],Zm[2*CDIMs];
  h_type JDms[2*CDIMs],ZDm[2*CDIMs];
  h_type mcnr,Norm;
  h_type xp,xs;

  if(kind<1 || kind>2)
     {std::cout <<"MathuZn: The kind has not been specified correctly:STOP"; exit(0);}

  xp=0.5*h*exp(u);
  xs=0.5*h*exp(-u);
  Coefficients(typ,n,h,Bm,CDIM,mcnr,Norm);

  BesselJ(xs,2*CDIM,Jms,JDms);
  if(kind==1) BesselJ(xp,2*CDIM,Zm,ZDm);
  if(kind==2) BesselY(xp,2*CDIM,Zm,ZDm);

   return MathuZDn(typ,n,h,u,Bm,Jms,JDms,Zm,ZDm,CDIM);

}
h_type MathuSn(char typ,int n, h_type h, double v, short kind)
{
/*
   inputs:
   typ ='e' or 'o'
   n  = order  
   h  = parameter
   v  = argument in radians
   kind = 1 or 2
  
   output:
   The function returns the value of the required circumferential function


*/
  int CDIM=CDIMs-1;
  h_type Bm[2*CDIMs],Dm[2*CDIMs];
  h_type Vr,mcnr,Norm1,Norm2;

  if(kind<1 || kind>2)
     {std::cout <<"MathuSn: The kind has not been specified correctly:STOP"; exit(0);}

  Coefficients(typ,n,h,Bm,CDIM,mcnr,Norm1);

  if(kind==1)  Vr=MathuSn(typ,n,v,Bm,CDIM);
  if(kind==2)
    {
     CoefficientSec(typ,n,h,mcnr,Bm,CDIM,Dm,Norm2);
     Vr=MathuFn(typ,n,v,Dm,Norm2,Bm,CDIM);
    }

   return Vr;

}
h_type MathuSDn(char typ,int n, h_type h, double v, short kind)
{
  int CDIM=CDIMs-1;
  h_type Bm[2*CDIMs],Dm[2*CDIMs];
  h_type Vr,mcnr,Norm1,Norm2;

  if(kind<1 || kind>2)
     {std::cout <<"MathuSn: The kind has not been specified correctly:STOP"; exit(0);}

  Coefficients(typ,n,h,Bm,CDIM,mcnr,Norm1);

  if(kind==1)  Vr=MathuSDn(typ,n,v,Bm,CDIM);
  if(kind==2)
    {
     CoefficientSec(typ,n,h,mcnr,Bm,CDIM,Dm,Norm2);
     Vr=MathuFDn(typ,n,v,Dm,Norm2,Bm,CDIM);
    }

   return Vr;

}
/*==============================================================*/
/*----------------- General way of computing Zen and Zon ------------*/
/*
   This Function takes the coefficients Bm[], the Bessel functions
   Jms[]=Jm(0.5h exp(-u))  Zm[]=Zm(0.5h exp(u)) all of size CDIM
   Z={Jn,Yn,Hn}
   This is a general Mathieu function which uses the best convergent series
   for Mathieu functions. This function is suitable for computing;
   first kind, second kind, third kind and fourth kind radial Mathieu functions
   of course Zm will have to changed as appropriate
   for Jen or Jon, Zm[]=Jm[]
       Yen or Yon  Zm[]=Ym[]
       He(1)n      Zm[]=H(1)m[]  etc
  Note for Hen some variables need to be changed to complex.

  Make sure when typ='e' to pass the Coefficients Bem[]
  and when  typ='o' pass the Coefficients Bom[]
*/
h_type  MathuZn(char typ, int n, h_type Bm[], h_type Jms[], h_type Zm[], int CDIM)
{
 int  m,ni,p,s;
 h_type Sum,MathZ,Tmp;

    p=n%2;
    ni=(n-p)/2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"MathuZn: Type has not been specified correctly:STOP";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"MathuZn: n can not be zero for odd functions :STOP";
	exit(0);
     }
    s=ni;
    if(s>=CDIM) s=CDIM-1;
    Sum=0;
    if ( typ=='e' )
    {
	    Tmp=1/(Bm[2*s+p]+1e-50);
	    if (s==0 && n!=1) Tmp=0.5/(Bm[2*s+p]+1e-50);

		for (m=s;m<CDIM-s-p;m++)
		    {Sum +=signf(ni-m)*Bm[2*m+p]*(Jms[m-s]*Zm[m+s+p]+Jms[m+s+p]*Zm[m-s]);}
		for (m=0;m<s;m++)
		    {Sum +=signf(ni-s)*Bm[2*m+p]*(Jms[s-m]*Zm[m+s+p]+Jms[m+s+p]*Zm[s-m]);}

     }

     if ( typ=='o' )
	{
	   Tmp=1/(Bm[2*s+p]+1e-50);
		for (m=s;m<CDIM-s-p;m++)
		    {Sum +=signf(ni-m)*Bm[2*m+p]*(Jms[m-s]*Zm[m+s+p]-Jms[m+s+p]*Zm[m-s]);}
		for (m=0;m<s;m++)
		    {Sum +=signf(ni-s)*Bm[2*m+p]*(Jms[s-m]*Zm[m+s+p]-Jms[m+s+p]*Zm[s-m]);}

     }
	  MathZ=pow(0.5*pi,0.5)*Sum*Tmp;
	   return MathZ;
}
/*-----------------------------------------------------------------------*/
/*- Differential of the Radial Second Kind Even and Odd Mathieu Function -*/
/*
  This Function takes the coefficients Bm[], the Bessel functions
  Jms[]=Jm(0.5h exp(-u)) and JDms[]. Zm[]=Zm(0.5h exp(u)) and ZDm[].
   where Zm={Jm,Ym,Hm}
  sums the series giving  Ze'n(h,cosh u) if typ='e'
  or Zo'n(h, cosh u) if typ='o'
  Note that differentiation is wrt u.
*/
h_type  MathuZDn(char typ, int n, h_type h, double u, h_type Bm[], h_type Jms[], h_type JDms[], h_type Zm[], h_type ZDm[], int CDIM)
{

 int  m,ni,p,s;
 h_type Sum,MathZD,Tmp;
 double un,up;

    un=exp(-u);
    up=exp(u);

    p=n%2;
    ni=(n-p)/2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MathuZDn: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n"<<"MathuZDn: n can not be zero for odd functions :STOP"<<"\n";
     exit(0);
     }

    s=ni;
    if(s>=CDIM) s=CDIM-1;
    Sum=0;
    if ( typ=='e' )
     {
	Tmp=1/(Bm[2*s+p]+1e-50);
	if (s==0 && n!=1) Tmp=0.5/(Bm[2*s+p]+1e-50);

	 for(m=s;m<CDIM-s-p;m++)
	  {Sum +=signf(ni-m)*Bm[2*m+p]*(up*ZDm[m+s+p]*Jms[m-s]-un*Zm[m+s+p]*JDms[m-s]+up*ZDm[m-s]*Jms[m+s+p]-un*Zm[m-s]*JDms[m+s+p]);}

	 for(m=0;m<s;m++)
	  {Sum +=signf(ni-s)*Bm[2*m+p]*(up*ZDm[m+s+p]*Jms[s-m]-un*Zm[m+s+p]*JDms[s-m]+up*ZDm[s-m]*Jms[m+s+p]-un*Zm[s-m]*JDms[m+s+p]);}


      }

    if ( typ=='o' )
     {
	Tmp=1/(Bm[2*s+p]+1e-50);

	 for(m=s;m<CDIM-s-p;m++)
	  {Sum +=signf(ni-m)*Bm[2*m+p]*(up*ZDm[m+s+p]*Jms[m-s]-un*Zm[m+s+p]*JDms[m-s]-up*ZDm[m-s]*Jms[m+s+p]+un*Zm[m-s]*JDms[m+s+p]);}

	 for(m=0;m<s;m++)
	  {Sum +=signf(ni-s)*Bm[2*m+p]*(up*ZDm[m+s+p]*Jms[s-m]-un*Zm[m+s+p]*JDms[s-m]-up*ZDm[s-m]*Jms[m+s+p]+un*Zm[s-m]*JDms[m+s+p]);}

      }
	MathZD=pow(0.5*pi,0.5)*Sum*0.5*h*Tmp;
     return MathZD;
}
/*=======================================================================*/
/*------- Circumferential First kind Even and Odd Mathieu  Function ------*/
h_type  MathuSn(char typ, int n, double v, h_type Bm[], int CDIM)
{
 int  m,p;
 h_type Sum,MathuS;

    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MathuSn: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n"<<"MathuSn: n can not be zero for odd functions :STOP"<<"\n";
     exit(0);
     }

    if ( typ=='e' )
       {
	  Sum=Bm[p]*cos(p*v);
	  for(m=1;m<CDIM;m++)
	     { Sum=Sum+Bm[2*m+p]*cos((2*m+p)*v);}

	  MathuS=Sum;
       }
    if ( typ=='o'  )
       {
	  Sum=Bm[p]*sin(p*v);
	  for(m=1;m<CDIM;m++)
	      {Sum=Sum+Bm[2*m+p]*sin((2*m+p)*v);}

	   MathuS=Sum;
       }

       return MathuS;

}
/*--------------------------------------------------------------------*/

/*--Dif. of Circumferential First kind Even and Odd Mathieu  Function --*/
h_type  MathuSDn(char typ, int n, double v, h_type Bm[], int CDIM)
{
 int  m,p;
 h_type Sum,MathuSD;

    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MathuSDn: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n"<<"MathuSDn: n can not be zero for odd functions :STOP"<<"\n";
     exit(0);
     }

    if (typ=='e')
       {
	  Sum=-Bm[p]*p*sin(p*v);
	  for(m=1;m<CDIM;m++)
	      {Sum=Sum-(2*m+p)*Bm[2*m+p]*sin((2*m+p)*v);}

	   MathuSD=Sum;
	}

     if ( typ=='o' )
	{
	  Sum=Bm[p]*p*cos(p*v);
	  for(m=1;m<CDIM;m++)
	      { Sum=Sum+(2*m+p)*Bm[2*m+p]*cos((2*m+p)*v);}

	  MathuSD=Sum;
	}

	return MathuSD;

}
/*--------------------------------------------------------------------*/
/*====================================================================*/
/*--- Circumferential Second Kind Mathieu Function (Non-periodic) ---*/
/*---------- Fen() and Fon()  ------*/
/*
   Bm is Bem or Bom 1st kind Mathieu coefficients
   Dm is Dem or Dom 2nd kind Mathieu coefficients
   Norm second kind normalization constant

*/
h_type  MathuFn(char typ, int n, double v, h_type Dm[], h_type Norm, h_type Bm[], int CDIM)
{
 int  m,p;
 h_type Sum,MathuF;

    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MathuFn: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n" <<"MathuFn: n can not be zero for odd functions :STOP"<<"\n";
     exit(0);
     }

    if ( typ=='e' )
      {
       Sum=Dm[p]*sin(p*v);
       for(m=1;m<CDIM;m++)
	   {Sum=Sum+Dm[2*m+p]*sin((2*m+p)*v);}
      }

    if ( typ=='o'  )
      {
       Sum=Dm[p]*cos(p*v);
       for (m=1;m<CDIM;m++)
	   {Sum=Sum+Dm[2*m+p]*cos((2*m+p)*v);}
       }

     MathuF=Norm*(v*MathuSn(typ,n,v,Bm,CDIM)+Sum);

    return MathuF;
}
/*----------------------------------------------------------------------*/
/*- Dif. of Circumferential Second Kind Mathieu Function (Non-periodic)--*/
/*------- Fe'n() and F'n() -------*/
h_type  MathuFDn(char typ, int n, double v, h_type Dm[], h_type Norm, h_type Bm[], int CDIM)
{
 int  m,p;
 h_type Sum,MathuFD;

    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MathuFDn: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n"<<"MathuFDn: n can not be zero for odd functions :STOP"<<"\n";
     exit(0);
     }


    if ( typ=='e' )
     {
      Sum=p*Dm[p]*cos(p*v);
      for(m=1;m<CDIM;m++)
       {Sum=Sum+Dm[2*m+p]*(2*m+p)*cos((2*m+p)*v);}
     }

    if ( typ=='o'  )
     { Sum=-p*Dm[p]*sin(p*v);
      for (m=1;m<CDIM;m++)
	{Sum=Sum-Dm[2*m+p]*(2*m+p)*sin((2*m+p)*v);}
    }

    MathuFD=Norm*(MathuSn(typ,n,v,Bm,CDIM)+v*MathuSDn(typ,n,v,Bm,CDIM)+Sum);

    return MathuFD;
}
/*--------------------------------------------------------------------*/
/*------------ Radial First Kind Even and Odd Mathieu Function ------*/
/*
  THESE SERIES ARE NOT CONVERGENT FOR ALL Z-PLANE
  BUT ARE LEFT HERE FOR TESTING PURPOSES.

 This function takes the coefficients Bm[] and the Bessel functions
 Jm(h cosh u)  of size 2*CDIM
 then sums the series giving Jen(h,cosh u) if typ='e'
 or  Jon(h, cosh u) if typ='o'
 Make sure when typ='e' to pass the Coefficients Bem[]
 and when  typ='o' pass the Coefficients Bom[]
*/

h_type  MathuJn(char typ, int n, double u, h_type Bm[], h_type Jm[], int CDIM)
{
 int  m,ni,p;
 h_type Sum,MathuJ;


    p=n%2;
    ni=(n-p)/2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"MathuJn: Type has not been specified correctly:STOP";
      exit(1);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"MathuJn: n can not be zero for odd functions :STOP";
     exit(2);
     }

    if ( typ=='e' )
     {
	    Sum=signf(ni)*Bm[p]*Jm[p];
	    for (m=1; m<CDIM; m++)
		{ Sum=Sum+signf(ni-m)*Bm[2*m+p]*Jm[2*m+p]; }

	    MathuJ=pow(pi/2.,.5)*Sum;
     }
    if ( typ=='o' )
      {
	    Sum=signf(ni)*p*Bm[p]*Jm[p];
	    for ( m=1; m<CDIM; m++)
		{ Sum=Sum+signf(ni-m)*Bm[2*m+p]*Jm[2*m+p]*(2*m+p); }

	    MathuJ=pow(pi/2.,.5)*Sum*tanh(u);
      }

       return MathuJ;
}
/*-----------------------------------------------------------------------*/
/*- Differential of the Radial First Kind Even and  Odd Mathieu Function-*/
/*
  This Function takes the coefficients Hm(), the Bessel functions Jm
  and the Dif  of Bessel functions Jm'(h cosh u) and n,h,u.
  sums the series giving  Je'n(h,cosh u) if typ='e'
  or Jo'n(h, cosh u) if typ='o'
  Note that differentiation is wrt u.
*/
h_type  MathuJDn(char typ, int n, h_type h, double u, h_type Bm[], h_type Jm[], h_type JDm[], int CDIM)
{
 int  m,ni,p;
 h_type Sum,MathuJD;

    p=n%2;
    ni=(n-p)/2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"MathuJDn: Type has not been specified correctly:STOP";
      exit(1);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"MathuJDn: n can not be zero for odd functions :STOP";
     exit(2);
     }

    if ( typ=='e' )
     {
	    Sum=signf(ni)*Bm[p]*JDm[p];
	    for ( m=1; m<CDIM; m++)
		{ Sum=Sum+signf(ni-m)*Bm[2*m+p]*JDm[2*m+p];}

	    MathuJD=pow(pi/2.,.5)*Sum*h*sinh(u);
     }
    if ( typ=='o' )
      {
	    Sum=signf(ni)*p*Bm[p]*(tanh(u)*sinh(u)*h*JDm[p]+Jm[p]/pow(cosh(u),2));
	    for ( m=1 ; m<CDIM; m++)
	     { Sum=Sum+signf(ni-m)*Bm[2*m+p]*(tanh(u)*sinh(u)*h*JDm[2*m+p]+Jm[2*m+p]/pow(cosh(u),2))*(2.*m+p);}

	     MathuJD=pow(pi/2.,.5)*Sum;
       }

       return MathuJD;

}
/*-----------------------------------------------------------------------*/
/*=======================================================================*/
/*------------ Radial Second Kind Even and Odd Mathieu Function ----------*/
/*
  This function takes the coefficients Hm() && the Bessel functions
  Jms(0.5h exp(-u)) and Ym(0.5h exp(u)).
  Then sums the series giving Yen(h,cosh u) if typ='e'
  or Yon(h,cosh u) if typ='o'.
*/
h_type  MathuYn(char typ, int n, h_type Bm[], h_type Jms[], h_type Ym[], int CDIM)
{
  int  m,ni,p;
  h_type Sum,MathuY;

    p=n%2;
    ni=(n-p)/2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"MathuYn: Type has not been specified correctly:STOP";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"MathuYn: n can not be zero for odd functions :STOP";
     exit(0);
     }

    if ( typ=='e' )
    { if ( p==0 )
	{   Sum=signf(ni)*Bm[0]*Ym[0]*Jms[0];
		for (m=1;m<CDIM;m++)
		    {Sum=Sum+signf(ni-m)*Bm[2*m]*Ym[m]*Jms[m];}
	}

     if ( p==1 )
	{     Sum=signf(ni)*Bm[1]*(Jms[0]*Ym[1]+Jms[1]*Ym[0]);
		for (m=1;m<CDIM-1;m++)
			{Sum=Sum+signf(ni-m)*Bm[2*m+1]*(Jms[m]*Ym[m+1]+Jms[m+1]*Ym[m]);}
	}

	  MathuY=pow(pi/2.,.5)*Sum/(Bm[p]+1e-30);
     }

     if ( typ=='o' )
	{
	   Sum=signf(ni)*Bm[p]*(Jms[1-p]*Ym[1]-Jms[1]*Ym[1-p]);
	   for (m=1;m<CDIM-1;m++)
		 {Sum=Sum+signf(ni-m)*Bm[2*m+p]*(Jms[m-1+p]*Ym[m+1]-Jms[m+1]*Ym[m-1+p]);}

	  MathuY=pow(pi/2.,.5)*Sum/(Bm[2-p]+1e-30);
	 }
	   return MathuY;
}
/*-----------------------------------------------------------------------*/
/*- Differential of the Radial Second Kind Even and Odd Mathieu Function -*/
/*
  This Function takes the coefficients Bm[], the Bessel functions
  Jms[]=Jm(0.5h exp(-u)) and JDms[]. Ym[]=Ym(0.5h exp(u)) and  YDm[].
  Sums the series giving  Ye'n(h,cosh u) if typ='e'
  or Yo'n(h, cosh u) if typ='o'
  Note that differentiation is wrt u.
*/
h_type MathuYDn(char typ, int n, h_type h, double u, h_type Bm[], h_type Jms[], h_type JDms[], h_type Ym[], h_type YDm[], int CDIM)
{
 int  m,ni,p;
 h_type Sum,MathuYD;
 double up,un;

 un=exp(-u);
 up=exp(u);

    p=n%2;
    ni=(n-p)/2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MathuYDn: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n"<<"MathuYDn: n can not be zero for odd functions :STOP"<<"\n";
     exit(0);
     }


    Sum=0;
    if ( typ=='e' )
     {
	if ( p==0 )
	{
	   for ( m=0;m<CDIM;m++)
		{Sum=Sum+signf(ni-m)*Bm[2*m]*(up*YDm[m]*Jms[m]-un*Ym[m]*JDms[m]);}
	}

	if ( p==1 )
	{
	   for (m=0;m<CDIM-1;m++)
	    {Sum=Sum+signf(ni-m)*Bm[2*m+1]*(up*YDm[m+1]*Jms[m]-un*Ym[m+1]*JDms[m]+up*YDm[m]*Jms[m+1]-un*Ym[m]*JDms[m+1]);}
	}

	MathuYD=pow(pi/2.,.5)*Sum*.5*h/(Bm[p]+1e-30);
      }

    if ( typ=='o' )
     {
	Sum=signf(ni)*Bm[p]*(up*YDm[1]*Jms[0]-un*JDms[0]*Ym[1]+un*JDms[1]*Ym[0]-up*Jms[1]*YDm[0]);
	for ( m=1 ; m<CDIM-1; m++)
	  { Sum=Sum+signf(ni-m)*Bm[2*m+p]*(up*YDm[m+1]*Jms[m-1+p]-un*Ym[m+1]*JDms[m-1+p]+un*JDms[m+1]*Ym[m-1+p]-up*Jms[m+1]*YDm[m-1+p]);}

	MathuYD=pow(pi/2.,.5)*Sum*.5*h/(Bm[2-p]+1e-30);
      }

     return MathuYD;
}
/*------------------------------------------------------------------------*/
/*====================================================================*/
/*---------- Integrals of Sen and Son  -------------------------------*/
/*
     This uses Romberg integration Routine
     The routine is general only the function call is chaged
*/
h_type ISn(char typ,int ni, double vi, double ti, double ui, h_type Bm[], int CDIM)
{
  int K,M,N,N1,N2,TMP;
  double a,b,ACC,DX,R,X;
  h_type T[2][11],INTG,Z;

  a=vi-ti;b=vi+ti;    // Limits: a lower limit, b upper limit
  INTG =MathuSn(typ,ni,a,Bm,CDIM)*pow(pow(sinh(ui),2)+pow(sin(a),2),.5);
  INTG+=MathuSn(typ,ni,b,Bm,CDIM)*pow(pow(sinh(ui),2)+pow(sin(b),2),.5);      //Function CALL
  T[0][0]=0.5*(b-a)*INTG;

  N1=0;N2=1;K=0;ACC=1;

  while((ACC>1e-6 && K<10) || K<4)    // Max NUMBER OF ITERATIONS
   {
    K++;
    N=(int)pow(2,K); DX=(b-a)/N;Z=0;

   for (R=1;R<=N;R+=2)
	{
	 X=a+(R*DX);
	 Z +=MathuSn(typ,ni,X,Bm,CDIM)*pow(pow(sinh(ui),2)+pow(sin(X),2),.5);      // Function CALL
	}

   T[N2][0]=.5*T[N1][0]+DX*Z;   // VALUE OF T(K,0)

   for ( M=1 ; M<=K; M++)
    { T[N2][M]=(pow(4,M)*T[N2][M-1]-T[N1][M-1])/(pow(4,M)-1);}

   ACC=abs(T[N2][K]-T[N1][K-1]);
   TMP=N1;N1=N2;N2=TMP; //SWITCH ARRAY
  }

  INTG=T[N1][K];
  if (ACC>0.05) std::cout <<"INTERGRATION DID NOT CONVERGE";


   return INTG;

}

/*========================================================================*/
/*--------------------------------------------------------------------*/
#undef h_type
SHAR_EOF
fi # end of overwriting check
if test -f 'mathur.h'
then
	echo shar: will not over-write existing file "'mathur.h'"
else
cat << "SHAR_EOF" > 'mathur.h'

   double MathuZn(char typ,int n, double h, double u, short kind);
   double MathuZDn(char typ,int n, double h, double u, short kind);

   double MathuSn(char typ,int n, double h, double v, short kind);
   double MathuSDn(char typ,int n, double h, double v, short kind);

   double  MathuJn(char typ, int n, double u, double Bm[], double Jm[], int CDIM);
   double  MathuJDn(char typ, int n, double h, double u, double Bm[], double Jm[], double JDm[], int CDIM);
   double  MathuYn(char typ, int n, double Bm[], double Jms[], double Ym[], int CDIM);
   double  MathuYDn(char typ, int n, double h, double u, double Bm[], double Jms[], double JDms[], double Ym[], double YDm[], int CDIM);

   double  MathuZn(char typ, int n, double Bm[], double Jms[], double Zm[], int CDIM);
   double  MathuZDn(char typ, int n, double h, double u, double Bm[], double Jms[], double JDms[], double Zm[], double ZDm[], int CDIM);

   double  MathuSn(char typ, int n, double v, double Bm[], int CDIM);
   double  MathuSDn(char typ, int n, double v, double Bm[], int CDIM);
   double  ISn(char typ,int ni, double vi, double ti, double ui, double Bm[], int CDIM);
   double  MathuFn(char typ, int n, double v, double Gm[], double Norm, double Bm[], int CDIM);
   double  MathuFDn(char typ, int n, double v, double Gm[], double Norm, double Bm[], int CDIM);

   double MmathuSn(char typ,int n, double h, double v);
   double  MmathuSn(char typ, int n, double v, double Bem[], double Bom[], int CDIM);

SHAR_EOF
fi # end of overwriting check
if test -f 'mbsslr.cpp'
then
	echo shar: will not over-write existing file "'mbsslr.cpp'"
else
cat << "SHAR_EOF" > 'mbsslr.cpp'
#include <cmath>
#include <iostream>
#include <cstdlib>

/*
   Implementation of Modified Bessel functions based on numerical recipes.
   The routins give accuracy of over 14 decimal places for most cases.
   Accuracy can be checked using the Wronskians

*/

	  double BesselIn(int n, double x);
	  double BesselKn(int n, double x);
	  void bessik(double x, double xnu, double *Iv, double *Kv, double *IDv, double *KDv);

	  void BesselI(double z, int M, double bIn[]);
	  void BesselI(double z, int M, double bIn[], double IDn[]);
	  void BesselK(double z, int M, double bKn[]);
	  void BesselK(double z, int M, double bKn[], double KDn[]);

/*-----------------------------------------------------------------*/
double BesselIn(int n, double x)
{
	int j,jsum=0,m;
	double bi,bim,bip,tox,ans,Sum=0;
	double ACR=40.0,sign=1;
	double I0,K0,ID0,KD0,I1;

	if (x == 0.0)	{if (n==0) return 1; else return 0.0;}

	bessik(x,0,&I0,&K0,&ID0,&KD0);
	I1=ID0;
	if (n==0) return I0;
	if (n==1) return I1;

		tox=2.0/fabs(x);
		bip=ans=0.0;
		bi=1.0;
		m=(n+(int) sqrt(ACR*n));
		if(m%2==0) sign=-1;
		for (j=2*m;j>0;j--) {
			bim=bip+j*tox*bi;
			bip=bi;
			bi=bim;
			if (fabs(bi) > 1.0e10) {
				ans *= 1.0e-10;
				bi *= 1.0e-10;
				bip *= 1.0e-10;
				Sum *=1.0e-10;
			}
			if (jsum) {sign=-sign;Sum+=sign*bi;}
			jsum=!jsum;
			if (j == n) ans=bip;
		}
		Sum=2*Sum-bi;
		ans *= I0/bi;
		return x < 0.0 && (n & 1) ? -ans : ans;

}
/*---------------------------------------------------------------*/
double BesselKn(int n, double x)
{
	int j;
	double bk,bkm,bkp,tox;
	double I0,K0,ID0,KD0,K1;
	bessik(x,0,&I0,&K0,&ID0,&KD0);
	K1=-KD0;

	if(n==0) return K0;
	if(n==1) return K1;

	tox=2.0/x;
	bkm=K0;
	bk=K1;
	for (j=1;j<n;j++) {
		bkp=bkm+j*tox*bk;
		bkm=bk;
		bk=bkp;
	}
	return bk;
}
/*----------------------------------------------------------*/
/*		Arrays of Modified Bessel Functions         */
/*==========================================================*/

void BesselI(double z, int M, double bIn[])
{
	int j,m,N;
	double bi,bim,bip,toz,fact,sign;
	double ACR=72.0,Sum;
	double I0,K0,ID0,KD0;

	if (z == 0.0)
	     {
	      bIn[0]=1.0;
	      for(j=1;j<M;j++) bIn[j]=0.0;
	      return;
	     }


		toz=2.0/fabs(z);
		bip=0.0;
		for(m=0;m<M;m++)   bIn[m]=0.0;
		bi=1.0;Sum=0.0;
		N=(M+(int) sqrt(ACR*M));
		if(N%2) sign=1; else sign=-1;
		for (j=2*N;j>0;j--) {
			bim=bip+j*toz*bi;
			bip=bi;
			bi=bim;
			if (fabs(bi) > 1.0e10) {
				for(m=0;m<M;m++) bIn[m]=bIn[m]*1.0e-10;
				bi *= 1.0e-10;
				bip *= 1.0e-10;
				Sum *= 1.0e-10;
			}
			if (j<M) bIn[j]=bip;
			if (j%2==0) {if(j%4==0) Sum+=bip; else Sum -=bip;}
		}
		Sum=2*Sum+bi;
		I0=bi/Sum;
		if(fabs(z)>1) bessik(z,0,&I0,&K0,&ID0,&KD0);
		bIn[0]=I0;
		fact=bIn[0]/bi;
		sign=1.0;
		for(j=1;j<M;j++)
		{  if(z<0.0) sign=-sign;
		   bIn[j]=sign*bIn[j]*fact;
		}

}

void BesselI(double z, int M, double bIn[], double IDn[])
{
   int n;

   BesselI(z,M,bIn);
   IDn[0]=bIn[1];
   for(n=1;n<M-1;n++)
    { IDn[n]=0.5*(bIn[n-1]+bIn[n+1]);}
   if(z==0) IDn[M-1]=0; else IDn[M-1]=bIn[M-2]-(M-1.0)*bIn[M-1]/z;

	return;
}

/*---------------------------------------------------------*/
void BesselK(double z,int M,double bKn[])
{
	int j;
	double bk,bkm,bkp,toz;
	double I0,K0,ID0,KD0,K1;

	bessik(z,0,&I0,&K0,&ID0,&KD0);
	for (j=0;j<M;j++) bKn[j]=0;

	K1=-KD0;
	toz=2.0/z;
	bkm=K0;
	bk=K1;
	bKn[0]=bkm;
	bKn[1]=bk;
	for (j=1;j<M-1;j++) {
		bkp=bkm+j*toz*bk;
		bkm=bk;
		bk=bkp;
		bKn[j+1]=bkp;
//		if(bkp>1e250) j=M;
	}
	return;
}
void BesselK(double z, int M, double bKn[], double KDn[])
{
   int n;

   BesselK(z,M,bKn);
   KDn[0]=-bKn[1];
   for(n=1;n<M-1;n++)
    {KDn[n]=-0.5*(bKn[n-1]+bKn[n+1]);}

   KDn[M-1]=-bKn[M-2]-(M-1.0)*bKn[M-1]/z;

	return;
}
/*--------------------------------------------------------------*/
/*---------   Fractional Order Modified Bessel Functions -------*/
/*--------------------------------------------------------------*/
void bessik(double x, double xnu, double *Iv, double *Kv, double *IDv, double *KDv)
{
	void beschb(double x, double *gam1, double *gam2, double *gampl,
		double *gammi);
	void nrerror(char error_text[]);
	int i,l,nl;
	double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1,gam2,
		gammi,gampl,h,p,pimu,q,q1,q2,qnew,Ivl,Ivl1,Ivmu,IDv1,IDvl,
		Ivtemp,Kv1,Kvmu,Kvmup,Kvtemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2;

	double EPS=1.0e-16,FPMIN=1.0e-30,PI=3.141592653589793;
	int MAXIT=10000,XMIN=2;


	if (x <= 0.0 || xnu < 0.0) std::cout<<"\n bad arguments in bessik";
	nl=(int)(xnu+0.5);
	xmu=xnu-nl;
	xmu2=xmu*xmu;
	xi=1.0/x;
	xi2=2.0*xi;
	h=xnu*xi;
	if (h < FPMIN) h=FPMIN;
	b=xi2*xnu;
	d=0.0;
	c=h;
	for (i=1;i<=MAXIT;i++) {
		b += xi2;
		d=1.0/(b+d);
		c=b+1.0/c;
		del=c*d;
		h=del*h;
		if (fabs(del-1.0) < EPS) break;
	}
	if (i > MAXIT) {std::cout<<"x too large in bessik; try asymptotic expansion";exit(0);}
	Ivl=FPMIN;
	IDvl=h*Ivl;
	Ivl1=Ivl;
	IDv1=IDvl;
	fact=xnu*xi;
	for (l=nl;l>=1;l--) {
		Ivtemp=fact*Ivl+IDvl;
		fact -= xi;
		IDvl=fact*Ivtemp+Ivl;
		Ivl=Ivtemp;
	}
	f=IDvl/Ivl;
	if (x < XMIN) {
		x2=0.5*x;
		pimu=PI*xmu;
		fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
		d = -log(x2);
		e=xmu*d;
		fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
		beschb(xmu,&gam1,&gam2,&gampl,&gammi);
		ff=fact*(gam1*cosh(e)+gam2*fact2*d);
		sum=ff;
		e=exp(e);
		p=0.5*e/gampl;
		q=0.5/(e*gammi);
		c=1.0;
		d=x2*x2;
		sum1=p;
		for (i=1;i<=MAXIT;i++) {
			ff=(i*ff+p+q)/(i*i-xmu2);
			c *= (d/i);
			p /= (i-xmu);
			q /= (i+xmu);
			del=c*ff;
			sum += del;
			del1=c*(p-i*ff);
			sum1 += del1;
			if (fabs(del) < fabs(sum)*EPS) break;
		}
		if (i > MAXIT) {std::cout<<"bessk seIves failed to converge";exit(0);}
		Kvmu=sum;
		Kv1=sum1*xi2;
	} else {
		b=2.0*(1.0+x);
		d=1.0/b;
		h=delh=d;
		q1=0.0;
		q2=1.0;
		a1=0.25-xmu2;
		q=c=a1;
		a = -a1;
		s=1.0+q*delh;
		for (i=2;i<=MAXIT;i++) {
			a -= 2*(i-1);
			c = -a*c/i;
			qnew=(q1-b*q2)/a;
			q1=q2;
			q2=qnew;
			q += c*qnew;
			b += 2.0;
			d=1.0/(b+a*d);
			delh=(b*d-1.0)*delh;
			h += delh;
			dels=q*delh;
			s += dels;
			if (fabs(dels/s) < EPS) break;
		}
		if (i > MAXIT) {std::cout<<"bessik: failure to converge in cf2";exit(0);}
		h=a1*h;
		Kvmu=sqrt(PI/(2.0*x))*exp(-x)/s;
		Kv1=Kvmu*(xmu+x+0.5-h)*xi;
	}
	Kvmup=xmu*xi*Kvmu-Kv1;
	Ivmu=xi/(f*Kvmu-Kvmup);
	*Iv=(Ivmu*Ivl1)/Ivl;
	*IDv=(Ivmu*IDv1)/Ivl;
	for (i=1;i<=nl;i++) {
		Kvtemp=(xmu+i)*xi2*Kv1+Kvmu;
		Kvmu=Kv1;
		Kv1=Kvtemp;
	}
	*Kv=Kvmu;
	*KDv=xnu*xi*Kvmu-Kv1;
}

void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi)
{
	double chebev(double a, double b, double c[], int m, double x);
	double xx;
	int NUSE1=7,NUSE2=8;
	static double c1[] = {
		-1.142022680371172e0,6.516511267076e-3,
		3.08709017308e-4,-3.470626964e-6,6.943764e-9,
		3.6780e-11,-1.36e-13};
	static double c2[] = {
		1.843740587300906e0,-0.076852840844786e0,
		1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
		2.42310e-10,-1.70e-13,-1.0e-15};

	xx=8.0*x*x-1.0;
	*gam1=chebev(-1.0,1.0,c1,NUSE1,xx);
	*gam2=chebev(-1.0,1.0,c2,NUSE2,xx);
	*gampl= *gam2-x*(*gam1);
	*gammi= *gam2+x*(*gam1);
}
double chebev(double a, double b, double c[], int m, double x)
{
	double d=0.0,dd=0.0,sv,y,y2;
	int j;

	if ((x-a)*(x-b) > 0.0) {std::cout<<"x not in range in routine chebev";exit(0);}
	y2=2.0*(y=(2.0*x-a-b)/(b-a));
	for (j=m-1;j>=1;j--) {
		sv=d;
		d=y2*d-dd+c[j];
		dd=sv;
	}
	return y*d-dd+0.5*c[0];
}
SHAR_EOF
fi # end of overwriting check
if test -f 'mbsslr.h'
then
	echo shar: will not over-write existing file "'mbsslr.h'"
else
cat << "SHAR_EOF" > 'mbsslr.h'
	   double BesselIn(int n, double x);
	   double BesselI0(double x);
	   double BesselI1(double x);
	   double BesselKn(int n, double x);
	   double BesselK0(double x);
	   double BesselK1(double x);

	  void BesselI(double z, int M, double bIn[]);
	  void BesselI(double z, int M, double bIn[], double IDn[]);
	  void BesselK(double z, int M, double bKn[]);
	  void BesselK(double z, int M, double bKn[], double KDn[]);
SHAR_EOF
fi # end of overwriting check
if test -f 'mcnr.cpp'
then
	echo shar: will not over-write existing file "'mcnr.cpp'"
else
cat << "SHAR_EOF" > 'mcnr.cpp'
#include   <iostream>
#include   <cstdlib>
#include   <cmath>

/*
  The Author would be happy to help in any way he can
  in the use of these routins.

  Contact Address:
  CERI, KACST, P.O. Box 6086,
  Riyadh 11442,
  Saudi Arabia,
  Fax:966+1+4813764
  Email: alhargan@kacst.edu.sa

  Legal Matters:
  The program is here as is, no guarantees are given, stated or implied,
  you may use it at your own risk.

  last modified Oct 1998
*/
/*
  If you do not understand the routines do not worry, all you need is
  the use of the functions
  Coefficients(...) for first kind functions
  CoefficientSec(...) for second kind circumferential functions
  These take care of evaluating the MCNs of the other routines.

  For theoretical derivation see:
  The companion paper
  or
  F.A. Alhargan, "A complete method for the computations of Mathieu
  characteristic numbers of integer orders," SIAM Review, vol-38, No. 2,
  pp. 239-255, June 1996.

*/

const double pi=3.141592653589793;
inline double abs(double z)
     {return  fabs(z);}
inline double acosh(double z)
       {return log(z+sqrt(z*z-1));}
inline double asinh(double z)
       {return log(z+sqrt(z*z+1));}
inline double atanh(double z)
       {return 0.5*log((1+z)/(1-z));}

/*--------------------------------------------------------------*/
/*                     MCNR header file                         */
/*--------------------------------------------------------------*/


     double Estimatmcn(char typ,int n, double h);
     double Estimatmcnl(char typ,int n, double t);
     double nChain(char typ, int n, double t);


     double MCNRoot(double a, char typ, int n, double h, double accr, double &acc);
     double NewtonRaphsons(double a, char typ, int n, double h, double &acp);

     double Asympn(char typ, int n, double h);
     double Mcapn(int n, double h);
     double Approxn(char typ, int n, double  h);

   double Coefficients(char typ,int n, double h, double Bm[], int CDIM, double &mcnr, double &Norm);
   int CoefficientSec(char typ, int n, double h, double mcn, double Bm[], int CDIM, double Gm[], double &Norm);

/*========================================================================*/

/*--------------------------------------------------------------------*/
/*------------------------------------------------------------------*/
double MCNRoot(double a, char typ, int n, double h,double accr, double &acc)
{
     /*
       inputs:
       a    = the estimated MCN
       typ  = type of MCN e or o
       n    = order of the MCN
       h    = the parameter
       accr = the required accuracy


       ouput:
       The function returns the improved MCN with accuracy acc
       acc  = the obtained accuracy


     */

    int i,aprxt;
    double df,tst,c,est;
    double MCNimpro,t1,act;

    act=accr;
    aprxt=0;
    MCNimpro=1.0*n*n;

   /* if h is sufficiently small then approx is accurate enough */
    if ( fabs(MCNimpro-a)<1e-16 && fabs(h)<0.25*n ) aprxt=1;


    if(aprxt)
      {
	 MCNimpro=a;
	 NewtonRaphsons(MCNimpro,typ,n,h,acc);
      }


     if(!aprxt)
     {
       c=1;MCNimpro=a;i=1;t1=0;
       t1=NewtonRaphsons(MCNimpro,typ,n,h,c);MCNimpro=t1;est=c;
       while (c>act && i<60 && acc!=0.0)  /* Start Newton iteration */
	{
	 t1=NewtonRaphsons(MCNimpro,typ,n,h,c);        //   c=abs(f(x)/f'(x))
	 acc=fabs(MCNimpro-t1);
	 MCNimpro=t1;
	 i=i+1;
	}
       if (c>1e-6 && c>act) std::cout <<"MCNRoot: Newton Raphsons has not converged Acc=" <<c;
       df=fabs(a-MCNimpro);
       tst=df/(fabs(a)+1);
       if (tst>0.15)  /* Assuming that the guess value is close to root */
	{
	 if (fabs(h)>fabs(n/5.0))
	 {
	   std::cout <<"\n";
	   std::cout << "MCNRoot:";
	   std::cout << "Warning  mcn not accurate enough, Type= " <<typ <<", n=" << n <<", h=" << h <<"\n";
	   std::cout << "Estimated Accrcy=" << est << ", Improved Accrcy=" << c <<"\n";
	   std::cout << "Estimated MCN ="<<a <<", Improved root="<< MCNimpro <<"\n";
	   std::cout << "Test=" <<tst <<"\n";
	   std::cout << "Estimated was taken to be the correct MCN as the Improved is in error";
	   std::cout <<"\n";
	 }
	 MCNimpro=a;
	}
       acc=c;
     }


     return MCNimpro;
}
/*-----------------------------------------------------------------*/
/*----------------------------------------------------------------*/
/*-------------------------------------------------------------------*/
/* This algorithem is from  J. Wimp p83 */

double NewtonRaphsons(double a, char typ, int n,double h, double &acp)
{
    /*
      This algorithem is from  J. Wimp p83
      also see equations (47)-(52).
      This computes the a-f(a)/f'(a) for the MCNs
      inputs:
       a=estimated MCN, typ='e' or 'o'
       n=order, h=parameter
      outputs:
       improved MCN
       accuracy of the computation acp
    */
    int p,ex=0;
    double  r,M=1.2*n+15;
    double  d,ih2,h2,Y0,Y1,Yr,Yr1,Yr2,Yd0,Yd1,Ydr,Ydr1,Ydr2,f,fd;
    double  nr,dm,mr;

    M=int(M);
    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"NewtonRaphsons: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n" <<"NewtonRaphsons: n can not be zero for odd functions :STOP" <<"\n";
       exit(0);
     }

    if(h==0) {acp=0; return 1.0*n*n;}  // when h is zero, MCN=n*n

    d=a;
    if ( p==0 && typ=='o' ) ex=1;
    h2=0.5*h;h2=h2*h2;
    ih2=1.0/h2;
    Ydr2=0;Ydr1=0;Yr2=0;Yr1=1;
    for ( r=M ; r>=ex; r--)
     {
       mr=(2.0*r+2.0+p);
       mr=mr*mr;
       dm=(d-mr);
       dm=dm/h2;
       Yr=dm*Yr1-Yr2;
       Ydr=Yr1*ih2+dm*Ydr1-Ydr2;
       Yr2=Yr1;Yr1=Yr;
       Ydr2=Ydr1;Ydr1=Ydr;
       nr=fabs(Yr);
       if(nr>1e4)  /* Normalize to avoid overflow */
	   {
	     nr=1/nr;
	     Yr1=Yr1*nr;
	     Yr2=Yr2*nr;
	     Ydr1=Ydr1*nr;
	     Ydr2=Ydr2*nr;
	   }

     }

    Y0=Yr1;Y1=Yr2;
    Yd0=Ydr1;Yd1=Ydr2;

     if ( typ=='e'  )
      {
       if ( p==0 ) {f= 0.5*d*Y0-h2*Y1; fd=0.5*d*Yd0-h2*Yd1+0.5*Y0;}
       if ( p==1 ) {f= (d-1-h2)*Y0-h2*Y1; fd=(d-1-h2)*Yd0-h2*Yd1+Y0;}
      }
     else
      {
       if ( p==1 ) {f=(d-1+h2)*Y0-h2*Y1; fd=(d-1+h2)*Yd0-h2*Yd1+Y0;}
       if ( p==0 ) {f=(d-4)*Y0-h2*Y1; fd=(d-4)*Yd0-h2*Yd1+Y0;}
      }

   fd=fd+1e-90;  /* in case fd=zero so it does not crash */
   acp=fabs(f);
   r=fabs(f/fd);
   if (acp>r) acp=r;     //This value is returned f(x)/f'(x)

   return d-f/fd;
}
/*=======================================================================*/
/*-----------------------------------------------------------------------*/
double Mcapn(int n, double h)
{
   double R,n2,hp2,hp4,hp8,hp12,m3,m5;
    /*
       McLachlan (p17) Approximation for |h|~<1.4n for n>3
       Approximation valid for both an and bn.
       Note for higher values of n (n>100) the limit is |h|~<1.1n
       and may be less as n becomes large.

       Computes McLachlan estimate of MCN for n>3
       inputs:
       n=order, h=parameter
       output:
       R=estimated MCN


    */
   if(h==0) return 0;
   hp2=0.25*h*h;
   hp4=hp2*hp2;
   hp8=hp4*hp4;
   hp12=hp8*hp4;
   n2=n*n;
   m3=(n2-1)*(n2-1)*(n2-1);
   m5=m3*(n2-1)*(n2-1);
   R=n2+hp4/(2.0*(n2-1))+(5*n+7)*hp8/(32.0*m3*(n2-4));
   if(h>0.2*n) R=R+(9.0*n2*n2+58.0*n2+29)*hp12/(64.0*m5*(n2-4.0)*(n2-9));

   return R;

}

/*------------------------------------------------------------------*/
double Asympn(char typ, int n, double h)
{
/*
  These functions give the asymptotic values for a(n) and b(n)
  in the region  2n<=~ h <infinity  for a(n)   n>0, n=0  3<=h<infinity
  in the region  2n<=~ h <infinity for b(n)  n>1, n=1  3<=h<infinity
  See McLachlan p. 232
  inputs:
    typ='e' or 'o', n=order, h=parameter
  output:
    R=estimated MCN

*/

    double m,m2,m3,m4,m6;
    double tw10,tw14,tw16,tw20;
    double ih,ih2,ih3,ih4,ih5;
    double R;

    if ( typ=='e' ) m=2*n+1; else m=2*n-1;

    m2=m*m;m3=m2*m;m4=m*m3;m6=m4*m2;
    ih=1.0/h;
    ih2=ih*ih;
    ih3=ih2*ih;ih4=ih3*ih;ih5=ih4*ih;
    tw10=1024.0;tw10=1.0/tw10; // pow(2,-10)
    tw14=tw10*0.0625;
    tw16=tw14*0.25;
    tw20=tw16*0.0625;

    R= m*h - 0.5*h*h-0.125*(m2+1)-0.015625*(m3+3*m)*ih;
    R=R-(tw10)*(5*m4+34*m2+9)*(ih2);
    R=R-(tw14)*m*(33*m4+410*m2+405)*(ih3)-(tw16)*(63*m6+1260*m4+2943*m2+486)*(ih4);
    R=R-(tw20)*m*(527*m6+15617*m4+69001*m2+41607)*(ih5);

     return R;

}
/*--------------------------------------------------------------------*/

double Approxn(char typ, int n, double  h)
{
/*
   This function gives the approximate value for a(n) in the region
   of 0<h<0.8n for n>3, for n<3 the region is 0<h<4.
   Also it gives the approximate value for b(n) in the region
   of 0<h<0.8n for n>3, for n<3 the region is 0<h<4
   These are new approximations see
     F.A. Alhargan, "A complete method for the computations of Mathieu
     characteristic numbers of integer orders," SIAM Review, vol-38, No. 2,
     pp. 239-255, June 1996.
   They give good results for n<=3. However, for n>3 McLachlan's Appr.
   give better results and greater region for h.

  inputs:
    typ='e' or 'o', n=order, h=parameter
  output:
    R=approximate MCN

*/
   double  Q,R,T,S,c2,c1,c0,t1,M,th,w;

   if (h==0)  return n*n;
   if ( (typ=='e') )
     {
     if ( n==0 )
	   {
	     T=h*h*0.25;
	     R= 2-sqrt(4.0+2.0*T*T);  // equ (20b)
	     return R;
	   }

     if ( n==1)
	   {
	     T=h*h*0.25;
	     S=5.0*T*T-16.0*T+64.0;
	     Q=sqrt(S);
	     R=5.0+0.5*(T-Q);  // equ (21a)
	     R=R;
	     return  R;
	   }
     if ( n==2 )   // equ (20d)
	{
	 T=h*h*0.25;
	 c0=20.0*T*T;
	 c1=-48.0-3.0*T*T;
	 c2=-8;
       }

     if ( n==3 )   // equ (21b)
	{
	 T=h*h*0.25;
	 c0=T*T*(T+8.0);
	 c1=16.*T-128.0-2.*T*T;
	 c2=-T-8;
	}

     }

   if ( typ=='o' )
    {
     if ( n==1)
	  {
	    T=h*h*0.25;
	    R=5-0.5*(T+sqrt(64.0+16.0*T+5.0*T*T));  // equ (21c)
	    return R;
	  }

     if ( n==2 )
	  {
	    T=h*h*0.25;
	    R=10-sqrt(36+T*T);  // equ (21d)
	    return R;
	  }

     if ( n==3 )  // equ (21e)
      {
	T=h*h*0.25;
	c0=T*T*(8-T);
	c1=-16.0*T-128-2.0*T*T;
	c2=T-8;
      }

     }

     /* This approximation developed for general n>3 in similar manner
	as the above approximation, but was not mentioned in the paper
	as it is not accurate enough, development shown in my PhD thesis */
     if ( n>3 )
       {
	T=h*h*0.25;
	c0=8.0*T*T;
	c1=(16-2*T*T-16.0*n*n);
	c2=-8;
       }

     /*  Solution of the cubic x^3+c2 x^2+c1 x+c0=0  */
     Q=(3*c1-c2*c2)/9 ;R=(9*c2*c1-27*c0-2*c2*c2*c2)/54;
     w=Q*Q*Q+R*R;
     if ( w>=0 )
      {
       t1=R+sqrt(w);
       S=abs(t1)/t1*pow(abs(t1),1/3.);
       t1=R-sqrt(w);
       T=abs(t1)/t1*pow(abs(t1),1/3.);
      }
     else
      {
       M=pow(R*R+fabs(w),1/6.);
       th=atan(sqrt(abs(w))/R)/3.0+4.0*pi/3.0;
       S=M*cos(th);T=M*cos(th);
      }
     R=(3*S+3*T-c2)/3.0;

     if ( R<0 && h>0.2*n )
       {
	 R=Asympn(typ,n,h);  /* in this case return asymptotic value */
       }
     else
       R=n*n+fabs(R);

   return R;
}

/*====================================================================*/
/*---------------------------------------------------------------------*/
double Estimatmcn(char typ, int n, double h)
{
  /*
    This is the main routine for estimating MCNs
    it combines all other routines (asymptotic,approximate,Chaining) to
    obtain an estimate of an MCN for  the given inputs.

  inputs:
    typ='e' or 'o', n=order, h=parameter
  output:
    estimated MCN

  */

  double ps[4];   /* Contains positions for regions where
		 approximation is valid see Table 3 for n<4*/
  double MCNE,Yn,Xn,Av,t;


  if(h==0) return n*n;

  if (n<4)
   {
   if ( typ=='e' )
     { ps[0]=4;ps[1]=4;ps[2]=3.5;ps[3]=5;}
    else
     {ps[1]=4;ps[2]=4.5;ps[3]=5.25;}

    if ( fabs(h)<=ps[n]   )   return Approxn(typ,n,h);
    if ( fabs(h)>ps[n]    )    return Asympn(typ,n,h);
   }


   t=fabs(h)/(1.0*n);
   Xn=Mcapn(n,h);
   Yn=Asympn(typ,n,h);
   Av=0.5*(Xn+Yn);

  if(Yn>Xn && t>=2.0) return Yn;
  if(Xn>Yn && t<=1.0)  return Xn;

  if (n<19)
   {
    if (t<=1.4)  return Xn;
    if (t>1.4)   return Yn;
   }


  if(n>18 && n<70)        /* note even can be  to n<79 & odd  to n<122 */
   {
     if(Xn>Yn && t<1.4) return Xn;
     if(Xn>Yn && t>1.7) return Yn;
     return Av;
   }


   MCNE=nChain(typ,n,t);       /* if the above does not cover the rang,
				  then  chain in n-direction */

   return MCNE;

}
/*----------------------------------------------*/
double nChain(char typ, int n, double t)
{
   double rt,ac,co,cn,h;
   int r,ro,m,p;
   p=n%2;
   m=50+p;
   ro=m+10;
   rt=1.0*ro/m;
   co=rt*rt*Estimatmcnl(typ,m,t);
   cn=co;

  for(r=m+50;r<n;r+=50)  /* it works for steps up to 50 */
   {
     rt=1.0*r/ro;
     h=t*ro;
     cn=rt*rt*MCNRoot(co,typ,ro,h,0.001,ac);
     co=cn;
     ro=r;
   }
     rt=1.0*n/ro;
     h=t*ro;
     cn=rt*rt*MCNRoot(co,typ,ro,h,0.001,ac);


     return cn;

}
/*--------------------------------------------------------------------*/
double Estimatmcnl(char typ, int n, double t)
{
  double MCNE,Yn,Xn,Av,h;

  if(n<20) {std::cout <<"Estimatmcnl: n can not be less than 20, n="<<n; exit(0);}
  if(n>60) {std::cout <<"Estimatmcnl: n can not be greater than 60, n="<<n; exit(0);}

   h=t*n;

   Xn=Mcapn(n,h);
   Yn=Asympn(typ,n,h);
   Av=0.5*(Xn+Yn);

   if(Xn>Yn && t<1.6) MCNE=Xn;
   if(Xn>Yn && t>1.6) MCNE=Yn;

   if(t>1.7) MCNE=Yn;
   if(t<1.5) MCNE=Xn;

   if(t>=1.5 && t<=1.7) MCNE=Av;

   MCNE=MCNRoot(MCNE,typ,n,h,0.001,Av);

   return MCNE;

}
/*---------------------------------------------------------------------*/

/*=====================================================================*/
/*============== Coeffecients For Mathieu Functions  ==================*/
/*
  This calculates the coefficients acoording to the standard porcedure
  outlined by McLachlan, with some modifications.
  The MCN is estimated then improved using MCNRoot.
  Since coefficients are sensitive to the accuracy of the MCN
  The main causes of non-convergence are
  1) Inaccurrate estimation of MCN
  2) The MCN is very close to a singularity
  3) problem 2) is overcome using forward and backward computation
  The  printed warnings often not very serious,
       however results must be checked
  The coefficient are normalized according to Morse

  Input: typ,n,h
  Output: Cofficients and Normlization constant M and
  accuracy of computation is returned


*/

/*-------- Coeffecients For Even Mathieu functions Of odd Order  B(2m+1) ----*/

double Coefficients(char typ,int n, double h, double Bm[], int CDIM, double &mcnr, double &Norm)
{
/*
       inputs:
       typ  = type of the coefficiens 'e' or 'o'
       n    = order of  the function
       h    = the parameter
       CDIM = number of the coefficients array Bm[]


       ouput:
       An array of the computed coefficients  Bm[]
       The function returns the accuracy of the computation
       mcnr = Obtained Mathieu characteristic Number
       Norm = Normalization constant
    


*/
  int MD,m,p,ni,ex=0;
  double acc,aprtest;
  double *V,Nrm,v0,Ti,Vni,mcn,h2,ih2;


    p=n%2;
    h2=0.25*h*h;

    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"Coefficients: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     {
      // std::cout <<"\n" <<"Coefficients: n can not be zero for odd functions :STOP" <<"\n"; exit(0);
	  for(m=0;m<=MD;m++)  Bm[2*m+p]=0;
	  mcnr=0;
	  Norm=1e250;   // infinity
	 return 0;
     }

     V= new double [CDIM];
    if (typ=='o' && p==0)  {ex=1;Bm[0]=0;}    else ex=0;

    MD=CDIM-1;

    mcn=Estimatmcn(typ,n,h);      //Estimate MCN={an,bn}
    acc=1e-8;
    mcn=MCNRoot(mcn,typ,n,h,1e-16,acc); // Improve the Accuracy of MCN
    mcnr=mcn;

      if (abs(h)==0)
      {
	for(m=0;m<=MD;m++)
	 {Bm[2*m+p]=0;}
	if ( typ=='e' )
	  {
	     // make sure that n<2*CDIM+p;
	    if(n<2*MD+p) Bm[n]=1;
	    if ( n==0 ) Nrm=2*pi; else Nrm=pi;
	   }
	else
	  { Bm[n]=1./n;Nrm=pi*(Bm[n]*Bm[n]); }

       Norm=Nrm;
       delete[] V;
       return 1e-8;
     }

    aprtest=fabs(mcn-n*n);

    if(aprtest<1e-7 && abs(h)<0.4*(n+1))   // Use approximation
      {
	for(m=0;m<=MD;m++) {Bm[2*m+p]=0;}

	if ( typ=='e' )
	  {
	   Bm[n]=1;
	   if ( n==0 ) Nrm=2*pi; else Nrm=pi;
	  }
	else
	 {
	   Bm[n]=1.0/n;
	   Nrm=pi*Bm[n]*Bm[n];
	 }

       Norm=Nrm;
       delete[] V;
       return 1e-6;
     }

  ih2=1.0/h2;
  if (n==0)  /* the n=0 is a special case where forward is not needed */
     {
      V[MD]=-h2/((2*MD+2+p)*(2*MD+2+p)-mcn);
      for(m=MD;m>=2-p;m--)
	 {
	  Ti=1.0/((2.0*m+p)*(2.0*m+p));
	  V[m-1]=-Ti*h2/(1-mcn*Ti+Ti*h2*V[m]);
	 }
      V[0]=-.5*h2/(1-mcn/4+h2*V[1]/4);v0=mcn*ih2;
      acc=abs(V[0]/v0-1);
    }
  else
    {   /* Use Forward & Backward  */

       ni=(n-p)/2;    /* Work out the point where f&b computation should meet */
       if ( ni>=MD ) ni=MD;

		     /* initialized for forward computation */
       if ( typ=='e'  )
       {
	if ( p==0 ) {V[0]= mcn*ih2;V[1]=(mcn-4)*ih2-2/V[0];}
	if ( p==1 ) {V[0]=((mcn-1)*ih2-1);V[1]=(mcn-9)*ih2-1/V[0];}
       }

       if ( typ=='o'  )
       {
	if ( p==0 ) {V[0]=0;V[1]=(mcn-4)*ih2;}
	if ( p==1 ) {V[0]=(1+(mcn-1)*ih2);V[1]=(mcn-9)*ih2-1/V[0];}
       }

       for(m=2;m<=ni;m++)  /* Forward Computation */
	 {V[m]=(mcn-(2.0*m+p)*(2.0*m+p))*ih2-1.0/V[m-1];}

      Vni=V[ni];  /* Store last value of forward computation */

		 /* start backward computation */
      V[MD]=-h2/((2*MD+2+p)*(2*MD+2+p)-mcn);  /* The previous V assumed to be zero */
      for(m=MD;m>ni;m--)
	{
	 Ti=1.0/((2.0*m+p)*(2.0*m+p));
	 V[m-1]=-Ti*h2/(1-mcn*Ti+Ti*h2*V[m]);
	}

      acc=abs(V[ni]/Vni-1);     /*Check that forward & backward agree*/
    }

      Bm[p+2*ex]=1;
      for(m=0+ex;m<MD;m++)
       {Bm[2*m+2+p]=V[m]*Bm[2*m+p];}

  /*- Normalization is different for the even and the odd Coeff. ----*/

    if ( typ=='e'  )
     {
       Nrm=1;
       for(m=1;m<=MD;m++)
	 {Nrm=Nrm+(Bm[2*m+p]);}

       Bm[p]=(1.0/Nrm);
       if ( p==0 ) Nrm=2*Bm[0]*Bm[0]; else Nrm=Bm[1]*Bm[1];
     }

   if ( typ=='o'  )
    {
       if ( p==0 ) Nrm=2; else Nrm=1;
       for(m=2-p;m<=MD;m++)
	  {Nrm=Nrm+(2.0*m+p)*Bm[2*m+p];}

       Bm[2-p]=(1.0/Nrm);
       Nrm=(Bm[2-p])*(Bm[2-p]);
    }

    if (acc>1e-5){ std::cout <<"\n";
		 std::cout << "Coefficients: Warning Coeffecient  Have not converged" <<"\n";
		 std::cout << "type is " << typ;
		 std::cout << " n=" << n << " h=" <<h;
		 std::cout << " acc=" << acc <<"\n";
	       }
   /*-------------------------------*/

    for(m=1+ex;m<=MD;m++)
      {
	Bm[2*m+p]=Bm[2*m+p]*Bm[p+2*ex];
	Nrm=Nrm+Bm[2*m+p]*Bm[2*m+p];
      }

    Norm=pi*Nrm;   //Normalization Constant

    delete[] V;
    return acc;

}

/*--------------------------------------------------------------------*/

int CoefficientSec(char typ, int n, double h, double mcn, double Bm[], int CDIM, double Gm[], double &Norm)
{
/*
   Calculation of the Second Kind Mathieu Coefficients 
   Inputs:
   Bm(),mcn,n,h,typ
   Bm= Coefficients of the first kind Bem or Bom
   mcn= Mathieu Characteristic Number an or bn
   n=order h=paramtere , typ= 'e' for even typ='o' for odd

   Outputs:
   Gm()=Second Kind Mathieu Coefficients,  Dem() or Dom()
   Norm=Second Kind Normalization constant
   returns 1 if computation has converged or 0 if not
*/

  int i,s,p,r,sign=1,conv=1;
  double  *w,h2,ih2,th,Nrm,acc;

    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"CoefficientSec: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     {
      // std::cout <<"\n" <<"CoefficientSec: n can not be zero for odd functions :STOP" <<"\n";exit(0);
	   for(r=0;r<=CDIM-1;r++)
	      Gm[2*r+p]=0;
       Norm=1e250;  // infinity
     }

    w= new double [2*CDIM+3];

   s=CDIM-1;h2=0.25*h*h;ih2=1.0/h2;

   if(typ=='o') sign=-1;
   w[2*s+p]=0;
   w[2*s+p-2]=-2.0*(2.0*s+p)*sign*Bm[2*s+p]*ih2;

   for(r=s-1;r>=1;r--)
     w[2*r+p-2]=((mcn-(2.0*r+p)*(2.0*r+p))*w[2*r+p]-2.0*(2.0*r+p)*sign*Bm[2*r+p])*ih2-w[2*r+2+p];

   if(typ=='e')
    {
     if ( p==0 )  th=((mcn-4.0)*w[2]*ih2-w[4])/(2.0*Bm[0])-2.0*mcn/(h2*h2);
     if ( p==1 )  th=((mcn-1+h2)*w[1]*ih2-w[3])/(2.0*Bm[1])-ih2;
    }

   if ( typ=='o' )
    {
     if ( p==0 )  {w[0]=(4.0*Bm[2]+(mcn-4.0)*w[2])/(2*h2)-0.5*w[4];th=(w[2]-mcn*w[0]*ih2)/Bm[2];}
     if ( p==1 )  th=(w[3]-(mcn-1-h2)*w[1]*ih2)/(2.0*Bm[1])-ih2;
    }

   for(r=0;r<=s;r++)
      Gm[2*r+p]=w[2*r+p]-th*Bm[2*r+p];

  if(typ=='e')
    {
     Gm[0]=0;
     Nrm=1;
     for(i=0;i<=s;i++)
      {Nrm=Nrm+(2.0*i+p)*Gm[2*i+p];}

     Nrm=1.0/Nrm;
    }

  if ( typ=='o' )
    {
     Nrm=Gm[p];
     for (i=1; i<=s; i++)
      { Nrm=Nrm+Gm[2*i+p];}

     Nrm=1/Nrm;
    }
    r=4;
    acc=(mcn-(2.0*r+p)*(2.0*r+p))*Gm[2*r+p]-h2*h2*(Gm[2*r+2+p]+Gm[2*r-2+p]);
    if(typ=='e') acc=acc-2.0*(2.0*r+p)*Bm[2*r+p];
	else acc=fabs(acc+2.0*(2.0*r+p)*Bm[2*r+p]);

    Norm=Nrm;
    if (acc>1e-5){ std::cout <<"\n";
		 std::cout << "CoefficientSec: Warning Coeffecient  Have not converged" <<"\n";
		 std::cout << "type is " << typ;
		 std::cout << " n=" << n << " h=" <<h;
		 std::cout << " acc=" << acc <<"\n";
		 conv=0;
	       }

    delete[] w;
    return conv;

}
/*==================================================================*/

/*====================================================================*/

SHAR_EOF
fi # end of overwriting check
if test -f 'mcnr.h'
then
	echo shar: will not over-write existing file "'mcnr.h'"
else
cat << "SHAR_EOF" > 'mcnr.h'

     double McLachlan(char typ,int n, double h);

     double Estimatmcn(char typ,int n, double h);
     double Estimatmcnl(char typ,int n, double t);
     double nChain(char typ, int n, double t);


     double MCNRoot(double a, char typ, int n, double h, double accr, double &acc);
     double NewtonRaphsons(double a, char typ, int n, double h, double &acp);

     double Asympn(char typ, int n, double h);
     double Mcapn(int n, double h);
     double Approxn(char typ, int n, double  h);

     double Coefficients(char typ,int n, double h, double Bm[], int CDIM, double &mcnr, double &Norm);
     int CoefficientSec(char typ, int n, double h, double mcn, double Bm[], int CDIM, double Gm[], double &Norm);
SHAR_EOF
fi # end of overwriting check
if test -f 'mmathur.cpp'
then
	echo shar: will not over-write existing file "'mmathur.cpp'"
else
cat << "SHAR_EOF" > 'mmathur.cpp'
#include <iostream>
#include <cstdlib>
#include <cmath>

/*========================================================================*/
/*  		         Modified Mathieu Functions                       */
/*========================================================================*/

/*
  The Author would be happy to help in any way he can
  in the use of these routins.

  CERI, KACST, P.O. Box 6086, Riydh, Saudi Arabia,
  Fax:966+1+4813764
  Email: alhargan@kacst.edu.sa


*/



/*
  Modified Mathieu functions of real argument
  Qn=(Qe,Qo,Ee,Eo) and MZ=(Ie,Io,Ke,Ko)

  last modified  Oct. 1998

*/
#define CDIMs 60  /* Important: CIDMs must be greater than the order n, i.e. CDIMs=n+15 */

const double pi=3.141592653589793;

inline double abs(double z)
     {return  fabs(z);}
inline double acosh(double z)
       {return log(z+sqrt(z*z-1));}
inline double asinh(double z)
       {return log(z+sqrt(z*z+1));}
inline double atanh(double z)
       {return 0.5*log((1+z)/(1-z));}
inline int signf(int m)
{
  if(m%2==0) return 1;
   else return -1;
}

#define h_type double
	/*
	  For complex h replace double by complex, routins in this file
	  do not need modifing. However, you will not be able to run the
	  complex  part unless you have the complex routins for computing
	  Mathieu  Charactristic Numbers (MCNs) mcnc.cpp.
	  The one I have just developed  has limited range for h complex,
	  I am working in its extension.
	  Though it will not be  available  for sometime.
	*/

     /* routines required from file mcnr.cpp  */
   h_type Estimatmcn(char typ,int n, h_type h);
   h_type MCNRoot(h_type a, char typ, int n, h_type h, double accr, double &acc);
   double Coefficients(char typ,int n, h_type h, h_type Bm[], int CDIM, h_type &mcnr, h_type &Norm);
   int CoefficientSec(char typ, int n, h_type h, h_type mcn, h_type Bm[], int CDIM, h_type Dm[], h_type &Norm);

   /* routines required from file bsslr.cpp */
       double BesselJ(double z, int M, double bJn[]);
       double BesselJ(double z, int M, double bJn[], double JDn[]);

       double BesselY(double z, int M, double Yn[]);
       double BesselY(double z, int M, double Yn[], double YDn[]);
   /* routines required from file mbsslr.cpp */
	  void BesselI(double z, int M, double bIn[]);
	  void BesselI(double z, int M, double bIn[], double IDn[]);
	  void BesselK(double z, int M, double bKn[]);
	  void BesselK(double z, int M, double bKn[], double KDn[]);

   /* routines required from file mathur.cpp */
   h_type  MathuSn(char typ, int n, double v, h_type Bm[], int CDIM);
   h_type  MathuSDn(char typ, int n, double v, h_type Bm[], int CDIM);
   h_type  MathuFn(char typ, int n, double v, h_type Dm[], h_type Norm, h_type Bm[], int CDIM);
   h_type  MathuFDn(char typ, int n, double v, h_type Dm[], h_type Norm, h_type Bm[], int CDIM);

   /*   Routines available in this file  */
   h_type MathuMZn(char typ,int n, h_type h, double u, short kind);
   h_type MathuMZDn(char typ,int n, h_type h, double u, short kind);

   h_type MathuQn(char typ,int n, h_type h, double v, short kind);
   h_type MathuQDn(char typ,int n, h_type h, double v, short kind);

   h_type  MathuIn(char typ, int n, h_type Am[], h_type Ims[], h_type Zm[], int CDIM);
   h_type  MathuKn(char typ, int n, h_type Am[], h_type Kms[], h_type Zm[], int CDIM);
   h_type  MathuIDn(char typ, int n, h_type h, double u, h_type Am[], h_type Ims[], h_type IDms[], h_type Zm[], h_type ZDm[], int CDIM);
   h_type  MathuKDn(char typ, int n, h_type h, double u, h_type Am[], h_type Ims[], h_type IDms[], h_type Zm[], h_type ZDm[], int CDIM);

   h_type  MathuQn(char typ, int n, double v, h_type Am[], int CDIM);
   h_type  MathuQDn(char typ, int n, double v, h_type Am[], int CDIM);
   h_type  IQn(char typ,int ni, double vi, double ti, double ui, h_type Am[], int CDIM);
   h_type  MathuEn(char typ, int n, double v, h_type Cm[], h_type Norm, h_type Am[], int CDIM);
   h_type  MathuEDn(char typ, int n, double v, h_type Cm[], h_type Norm, h_type Am[], int CDIM);


   double MCoefficients(char typ,int n, h_type h, h_type Am[], int CDIM, h_type &mcnr, h_type &Norm);
   int MCoefficientSec(char typ, int n, h_type h, h_type mcn, h_type Am[], int CDIM, h_type Cm[], h_type &Norm);

/* ====================================================================*/
/* ================== Mathieu Functions ===============================*/
/*
  Important Convensions
1. CDIM is the number of terms summed in the series, which also determins
   the size of the arrays.
2. The value p is used to identify if n is even (p=0) or odd (p=1).
3. D in any function means differential e.g JDn(x)= d(Jn(x))/dx=J'n(x)
4. The series used are that which are convergent in all the z-plane.

5. The Coefficients of the Even functions is assumed be stored in
   one array Be() for both even and odd orders, similarly for odd functions Bo().
6. The Coefficients Be() and Bo() have to be calculated before calling the routine
   for Mathieu functions.
7. In the case of Jen and Jon the Bessel functions has to be calculated first
   Also for Yen and Yon.
8. Warning are not serious, but care needs to be taken as to validity of the
   results, in most cases the results are correct.
9. In the case of a serious error the programs haults with a messag

10.  Bm stands for Bem or Bom     1st kind coefficients
     Dm stands for Dem or Dom     2nd kind coefficients

11. All arrays are of size 2*CDIM, i.e Bm[0]...Bm[2*CDIM-1]
    some arrays are of size CDIM, e.g. Zm[0]...Zm[CDIM-1]

*/
/*------------------------------------------------------------------*/
/*---------------- Easy Functions -----------------------------------*/
/*
  When computing Mathieu functions, one usually would be requiring several
  different functions of the same order, in this case it is best to use
  the functions which take the coefficients and arrays of Bessel functions
  as inputs.
  On the other hand, one may want to compute a single function of one order
  in this case, 'The Easy Functions' coming shortly may be better and easily
  used.
*/
double MathuMZn(char typ,int n, double h, double u, short kind)
{
/*
   inputs:
   typ ='e' or 'o'
   n  = order
   h  = parameter
   u  = argument 
   kind = 1 or 2

   output:
   The function returns the value of the required modified radial function


*/

  int CDIM=CDIMs-1;
  double Am[2*CDIMs],Ims[2*CDIMs],Zm[2*CDIMs];
  double mcnr,Norm,Vr;
  double xp,xs;
  

  if(kind<1 || kind>2)
     {std::cout <<"MathuMZn: The kind has not been specified correctly:STOP"; exit(0);}


  xp=0.5*h*exp(u);
  xs=0.5*h*exp(-u);

  MCoefficients(typ,n,h,Am,CDIM,mcnr,Norm);

  BesselI(xs,2*CDIM,Ims);
  if(kind==1)
    { BesselI(xp,2*CDIM,Zm);
      Vr=MathuIn(typ,n,Am,Ims,Zm,CDIM);
    }
  if(kind==2)
     {  BesselK(xp,2*CDIM,Zm);
	Vr=MathuKn(typ,n,Am,Ims,Zm,CDIM);
     }

      return Vr;


}
double MathuMZDn(char typ,int n, double h, double u, short kind)
{
  int CDIM=CDIMs-1;
  double Am[2*CDIMs],Ims[2*CDIMs],Zm[2*CDIMs];
  double IDms[2*CDIMs],ZDm[2*CDIMs];
  double mcnr,Norm,Vr;
  double xp,xs;
  

  if(kind<1 || kind>2)
     {std::cout <<"MathuMZn: The kind has not been specified correctly:STOP"; exit(0);}


  xp=0.5*h*exp(u);
  xs=0.5*h*exp(-u);


  MCoefficients(typ,n,h,Am,CDIM,mcnr,Norm);

  BesselI(xs,2*CDIM,Ims,IDms);
  if(kind==1)
    { BesselI(xp,2*CDIM,Zm,ZDm);
       Vr=MathuIDn(typ,n,h,u,Am,Ims,IDms,Zm,ZDm,CDIM);
    }
  if(kind==2)
     {  BesselK(xp,2*CDIM,Zm,ZDm);
	Vr=MathuKDn(typ,n,h,u,Am,Ims,IDms,Zm,ZDm,CDIM);
     }

    return Vr;
}
/*------------------------------------------------------*/
double MathuQn(char typ,int n, double h, double v, short kind)
{
/*
   inputs:
   typ ='e' or 'o'
   n  = order
   h  = parameter
   v  = argument in radians
   kind = 1 or 2

   output:
   The function returns the value of the required modified circumferential function


*/

  int CDIM=CDIMs-1;
  h_type Am[2*CDIMs],Cm[2*CDIMs];
  h_type Vr,mcnr,Norm1,Norm2;


  MCoefficients(typ,n,h,Am,CDIM,mcnr,Norm1);

  if(kind==1) Vr=MathuQn(typ,n,v,Am,CDIM);
  if(kind==2)
  {
   MCoefficientSec(typ,n,h,mcnr,Am,CDIM,Cm,Norm2);
   Vr=MathuEn(typ,n,v,Cm,Norm2,Am,CDIM);
  }

   return Vr;

}
double MathuQDn(char typ,int n, double h, double v, short kind)
{
  int CDIM=CDIMs;
  h_type Am[2*CDIMs],Cm[2*CDIMs];
  h_type Vr,mcnr,Norm1,Norm2;


  MCoefficients(typ,n,h,Am,CDIM,mcnr,Norm1);

  if(kind==1) Vr=MathuQDn(typ,n,v,Am,CDIM);
  if(kind==2)
  {
   MCoefficientSec(typ,n,h,mcnr,Am,CDIM,Cm,Norm2);
   Vr=MathuEDn(typ,n,v,Cm,Norm2,Am,CDIM);
  }

   return Vr;

}
/*=======================================================================*/
/*====================================================================*/
/*-------- Integrals of Qen && Qon For Ports away from the center -----*/
/*
     This uses Romberg integration Routine
     The routine is general only the function call is chaged
*/
h_type IQn(char typ,int ni, double vi, double ti, double ui, h_type Am[], int CDIM)
{
  int K,M,N,N1,N2,TMP;
  double a,b,ACC,DX,R,X;
  double sui,sna,snb;
  h_type T[2][11],INTG,Z;

  a=vi-ti;b=vi+ti;    // Limits
  sui=sinh(ui);
  sui=sui*sui;
  sna=sin(a); sna=sna*sna;
  snb=sin(b);snb=snb*snb;
  INTG =MathuQn(typ,ni,a,Am,CDIM)*sqrt(sui+sna);
  INTG+=MathuQn(typ,ni,b,Am,CDIM)*sqrt(sui+snb);      //Function CALL
  T[0][0]=0.5*(b-a)*INTG;

  N1=0;N2=1;K=0;ACC=1;
  N=1;
  while((ACC>1e-6 && K<10) || K<4)    // Max NUMBER OF ITERATIONS
   {
    K++;
    N=2*N; DX=(b-a)/N;Z=0.0;

   for (R=1;R<=N;R+=2)
	{
	 X=a+(R*DX);
	 Z +=MathuQn(typ,ni,X,Am,CDIM)*pow(pow(sinh(ui),2)+pow(sin(X),2),.5);      // Function CALL
	}

   T[N2][0]=.5*T[N1][0]+DX*Z;   // VALUE OF T(K,0)

   for ( M=1 ; M<=K; M++)
    { T[N2][M]=(pow(4.0,M)*T[N2][M-1]-T[N1][M-1])/(pow(4.0,M)-1);}

   ACC=abs(T[N2][K]-T[N1][K-1]);
   TMP=N1;N1=N2;N2=TMP; //SWITCH ARRAY
  }

  INTG=T[N1][K];
  if (ACC>0.05) std::cout <<"INTERGRATION DID NOT CONVERGE";


   return INTG;

}
/*--------------------------------------------------------------------*/
/*====================================================================*/
/*---Modified Circumferential First kind Even && Odd Mathieu  Function ---*/
double  MathuQn(char typ, int n, double v, double Am[],  int CDIM)
{

 double MathuQ;
		/*
		  Make sure that the coefficients correspond correctly to
		  the type and the order.
		  i.e the coefficients passed are assumed to be
		  Modified Mathieu Coefficients;
		*/

      MathuQ=MathuSn(typ,n,v,Am,CDIM);

       return MathuQ;

}
/*--------------------------------------------------------------------*/
/*-Dif. of Modified Circumferential First kind Even && Odd Mathieu  Function-*/
double  MathuQDn(char typ, int n, double v, double Am[], int CDIM)
{

 double MathuQ;

	MathuQ=MathuSDn(typ,n,v,Am,CDIM);

       return MathuQ;

}
/*--- Circumferential Second Kind Mathieu Function (Non-periodic) ---*/
/*---------- Fen() && Fon()  ------*/
/*
   Bm is Bem or Bom 1st kind Mathieu coefficients
   Dm is Dem or Dom 2nd kind Mathieu coefficients
   Norm=second kind normalization constant
    not finished yet
*/
h_type  MathuEn(char typ, int n, double v, h_type Cm[], h_type Norm, h_type Am[], int CDIM)
{

 h_type MathuE;

	  /*   assuming that Cm,Am and Norm are the modified elements  */


     MathuE=MathuFn(typ,n,v,Cm,Norm,Am,CDIM);

    return MathuE;
}
/*----------------------------------------------------------------------*/
/*- Dif. of Circumferential Second Kind Mathieu Function (Non-periodic)--*/
/*------- Fe'n() and F'n() -------*/
h_type  MathuEDn(char typ, int n, double v, h_type Cm[], h_type Norm, h_type Am[], int CDIM)
{
 
 h_type MathuED;

	  /*   assuming that Cm,Am and Norm are the modified Coefficients  */

     MathuED=MathuFDn(typ,n,v,Cm,Norm,Am,CDIM);

    return MathuED;
}
/*--------------------------------------------------------------------*/
/*--------------------------------------------------------------------*/
/*------ Modified Radial Mahtieu function  -----------------*/
h_type  MathuIn(char typ, int n, h_type Am[], h_type Ims[], h_type Zm[], int CDIM)
{
 int  m,ni,p,s;
 h_type Sum,MathZ,Tmp;

	 /* assuming that  Am's are the Modified Mathieu Coefficients */

    p=n%2;
    ni=(n-p)/2;

    s=ni;
    if(s>=CDIM) s=CDIM-1;
    Sum=0;
    if (typ=='e')
    {
	    Tmp=1/(Am[2*s+p]+1e-50);
	    if (s==0 && n!=1) Tmp=0.5/(Am[2*s+p]+1e-50);

		for (m=s;m<CDIM-s-p;m++)
		    {Sum +=Am[2*m+p]*(Ims[m-s]*Zm[m+s+p]+Ims[m+s+p]*Zm[m-s]);}
		for (m=0;m<s;m++)
		    {Sum +=Am[2*m+p]*(Ims[s-m]*Zm[m+s+p]+Ims[m+s+p]*Zm[s-m]);}
     }

     if ( typ=='o' )
	{
	   Tmp=1/(Am[2*s+p]+1e-50);
		for (m=s;m<CDIM-s-p;m++)
		    {Sum +=Am[2*m+p]*(Ims[m-s]*Zm[m+s+p]-Ims[m+s+p]*Zm[m-s]);}
		for (m=0;m<s;m++)
		    {Sum +=Am[2*m+p]*(Ims[s-m]*Zm[m+s+p]-Ims[m+s+p]*Zm[s-m]);}
     }
	  MathZ=Sum*Tmp;
	   return MathZ;
}
/*------------------------------------------------------------------------------*/
h_type  MathuKn(char typ, int n, h_type Am[], h_type Ims[], h_type Zm[], int CDIM)
{
 int  m,ni,p,s;
 h_type Sum,MathZ,Tmp;
 double sign=1.0;

    p=n%2;
    ni=(n-p)/2;
    if(p==1) sign=-1.0;
    s=ni;
    if(s>=CDIM) s=CDIM-1;
    Sum=0;
    if ( typ=='e' )
    {
	    Tmp=1.0/(Am[2*s+p]+1e-50);
	    if (s==0 && n!=1) Tmp=0.5/(Am[2*s+p]+1e-50);

		for (m=s;m<CDIM-s-p;m++)
		    {Sum +=signf(m)*Am[2*m+p]*(Ims[m-s]*Zm[m+s+p]
						+sign*Ims[m+s+p]*Zm[m-s]);}
		for (m=0;m<s;m++)
		    {Sum +=signf(m)*Am[2*m+p]*(Ims[s-m]*Zm[m+s+p]
						+sign*Ims[m+s+p]*Zm[s-m]);}
     }

     if ( typ=='o' )
	{
	   Tmp=1/(Am[2*s+p]+1e-50);
		for (m=s;m<CDIM-s-p;m++)
		    {Sum +=signf(m)*Am[2*m+p]*(Ims[m-s]*Zm[m+s+p]
						-sign*Ims[m+s+p]*Zm[m-s]);}
		for (m=0;m<s;m++)
		    {Sum +=signf(m)*Am[2*m+p]*(Ims[s-m]*Zm[m+s+p]
						-sign*Ims[m+s+p]*Zm[s-m]);}
     }
	  MathZ=Sum*Tmp*signf(s);

   return MathZ;
}
/*--------------------------------------------------------------------*/

/*-----------------------------------------------------------------------*/
/*- Differential of the Radial Second Kind Even and Odd Mathieu Function -*/
/*
  This Function takes the coefficients Am[], the Bessel functions
  Ims[]=Im(0.5h exp(-u)) and IDms[]. Zm[]=Zm(0.5h exp(u)) and ZDm[].
   where Zm=Im
  sums the series giving  Ze'n(h,cosh u) if typ='e'
  or Zo'n(h, cosh u) if typ='o'
  Note that differentiation is wrt u.
*/
h_type  MathuIDn(char typ, int n, h_type h, double u, h_type Am[], h_type Ims[], h_type IDms[], h_type Zm[], h_type ZDm[], int CDIM)
{

 int  m,ni,p,s;
 h_type Sum,MathZD,Tmp;
 double un,up;

	/* assuming that Am's are the Modified Mathieu Coefficients */
    un=exp(-u);
    up=exp(u);

    p=n%2;
    ni=(n-p)/2;

    s=ni;
    if(s>=CDIM) s=CDIM-1;
    Sum=0;
    if ( typ=='e' )
     {
	Tmp=1/(Am[2*s+p]+1e-50);
	if (s==0 && n!=1) Tmp=0.5/(Am[2*s+p]+1e-50);

	 for(m=s;m<CDIM-s-p;m++)
	  {Sum +=Am[2*m+p]*(up*ZDm[m+s+p]*Ims[m-s]-un*Zm[m+s+p]*IDms[m-s]
			   +up*ZDm[m-s]*Ims[m+s+p]-un*Zm[m-s]*IDms[m+s+p]);}

	 for(m=0;m<s;m++)
	  {Sum +=Am[2*m+p]*(up*ZDm[m+s+p]*Ims[s-m]-un*Zm[m+s+p]*IDms[s-m]
			   +up*ZDm[s-m]*Ims[m+s+p]-un*Zm[s-m]*IDms[m+s+p]);}
      }

    if ( typ=='o' )
     {
	Tmp=1/(Am[2*s+p]+1e-50);

	 for(m=s;m<CDIM-s-p;m++)
	  {Sum +=Am[2*m+p]*(up*ZDm[m+s+p]*Ims[m-s]-un*Zm[m+s+p]*IDms[m-s]
			   -up*ZDm[m-s]*Ims[m+s+p]+un*Zm[m-s]*IDms[m+s+p]);}

	 for(m=0;m<s;m++)
	  {Sum +=Am[2*m+p]*(up*ZDm[m+s+p]*Ims[s-m]-un*Zm[m+s+p]*IDms[s-m]
			    -up*ZDm[s-m]*Ims[m+s+p]+un*Zm[s-m]*IDms[m+s+p]);}
      }

	MathZD=Sum*.5*h*Tmp;

     return MathZD;
}
/*------------------------------------------------------------------*/
h_type  MathuKDn(char typ, int n, h_type h, double u, h_type Am[], h_type Ims[], h_type IDms[], h_type Zm[], h_type ZDm[], int CDIM)
{

 int  m,ni,p,s;
 h_type Sum,MathZD,Tmp;
 double un,up,sign=1.0;

    un=exp(-u);
    up=exp(u);

    p=n%2;
    if(p==1) sign=-1.0;
    ni=(n-p)/2;
    s=ni;
    if(s>=CDIM) s=CDIM-1;
    Sum=0;
    if ( typ=='e' )
     {
	Tmp=1.0/(Am[2*s+p]+1e-50);
	if (s==0 && n!=1) Tmp=0.5/(Am[2*s+p]+1e-50);

	 for(m=s;m<CDIM-s-p;m++)
	  {Sum +=signf(m)*Am[2*m+p]*(up*ZDm[m+s+p]*Ims[m-s]-un*Zm[m+s+p]*IDms[m-s]
			       +sign*(up*ZDm[m-s]*Ims[m+s+p]-un*Zm[m-s]*IDms[m+s+p]));}

	 for(m=0;m<s;m++)
	  {Sum +=signf(m)*Am[2*m+p]*(up*ZDm[m+s+p]*Ims[s-m]-un*Zm[m+s+p]*IDms[s-m]
			       +sign*(up*ZDm[s-m]*Ims[m+s+p]-un*Zm[s-m]*IDms[m+s+p]));}

      }

    if ( typ=='o' )
     {
	Tmp=1.0/(Am[2*s+p]+1e-50);

	 for(m=s;m<CDIM-s-p;m++)
	  {Sum +=signf(m)*Am[2*m+p]*(up*ZDm[m+s+p]*Ims[m-s]-un*Zm[m+s+p]*IDms[m-s]
			       -sign*(up*ZDm[m-s]*Ims[m+s+p]-un*Zm[m-s]*IDms[m+s+p]));}

	 for(m=0;m<s;m++)
	  {Sum +=signf(m)*Am[2*m+p]*(up*ZDm[m+s+p]*Ims[s-m]-un*Zm[m+s+p]*IDms[s-m]
			       -sign*(up*ZDm[s-m]*Ims[m+s+p]-un*Zm[s-m]*IDms[m+s+p]));}

      }

      MathZD=Sum*.5*h*Tmp*signf(s);

     return MathZD;
}
/*=======================================================================*/
/*----------  Modified Mathieu Coefficients of the first kind ----------*/
double MCoefficients(char typ,int n, h_type h, h_type Am[], int CDIM, h_type &mcnr, h_type &Norm)
{
/*
   Calculation of the Modified First Kind  Mathieu Coefficients
   Inputs:
   n,h,typ
   n=order, h=paramtere , typ= 'e' for even typ='o' for odd

   Outputs:
   Am= First Kind Modified Mathieu  Coefficients  Aem or Aom
   mcn= Modified Mathieu Characteristic Number 
   Norm=First Kind Modified Normalization constant
   returns the accuracy of the computation acc


*/
  int p,m;
  char ntyp;
  double acc,Tmp;

  p=n%2;
  ntyp=typ;
  if(p==1) {if (typ=='e') ntyp='o'; else ntyp='e';}
  acc=Coefficients(ntyp,n,h,Am,CDIM,mcnr,Norm);

    if (typ=='e') Tmp=1.0/MathuSn(ntyp,n,0.5*pi,Am,CDIM);
    if (typ=='o') Tmp=1.0/MathuSDn(ntyp,n,0.5*pi,Am,CDIM);
    if (typ=='o' && p==1) Tmp=-Tmp;

    for(m=0;m<CDIM;m++)
      {
	   Am[2*m+p]=signf(m)*Am[2*m+p]*Tmp;
      }

    Norm=Norm*Tmp;

  return acc;
}
/*---------  Modified Second Kind Coeffecients --------------------*/
int MCoefficientSec(char typ, int n, h_type h, h_type mcn, h_type Am[], int CDIM, h_type Cm[], h_type &Norm)
{

	  /* The same routin for computing the standard 2nd Kind Coef
	     is used for computing the Modified 2nd Kind Coef
	     with small modifications
	     1. Bm replace by modified coeficients Am
	     2. h2 is replaced by -h2
	     3. for (p==1) mcn is interchanged between the even and the odd
	     Input:typ,n,h,mcn,Am[],CDIM
	     Output: Cm[],Norm
	   */

  int i,s,p,r,sign=1,conv=1;
  h_type  *w,h2,ih2,th,Nrm;

    w= new h_type [2*CDIM+3];
    p=n%2;
    if ( (typ!='e') && (typ!='o'))
    { std::cout <<"\n"<<"MCoefficientSec: Type has not been specified correctly: type=" << typ <<"\n";
      exit(0);
    }
    if ( n==0 && typ=='o' )
     { std::cout <<"\n" <<"MCoefficientSec: n can not be zero for odd functions :STOP" <<"\n";
     exit(0);
     }

     s=CDIM-1;
     h2=-0.25*h*h;  // it is negative for the modified
     ih2=1.0/h2;

   if(typ=='o') sign=-1;
   w[2*s+p]=0;
   w[2*s+p-2]=-2.0*(2.0*s+p)*sign*Am[2*s+p]*ih2;

   for(r=s-1;r>=1;r--)
     w[2*r+p-2]=((mcn-(2.0*r+p)*(2.0*r+p))*w[2*r+p]-2.0*(2.0*r+p)*sign*Am[2*r+p])*ih2-w[2*r+2+p];

   if(typ=='e')
    {
     if ( p==0 )  th=((mcn-4.0)*w[2]*ih2-w[4])/(2.0*Am[0])-2.0*mcn/(h2*h2);
     if ( p==1 )  th=((mcn-1.0+h2)*w[1]*ih2-w[3])/(2*Am[1])-ih2;
    }

   if ( typ=='o' )
    {
     if ( p==0 )  {w[0]=(4*Am[2]+(mcn-4)*w[2])/(2*h2)-w[4]/2;th=(w[2]-mcn*w[0]*ih2)/Am[2];}
     if ( p==1 )  th=(w[3]-(mcn-1-h2)*w[1]*ih2)/(2*Am[1])-ih2;
    }

   for(r=0;r<=s;r++)
      Cm[2*r+p]=w[2*r+p]-th*Am[2*r+p];

  if(typ=='e')
    {
     Cm[0]=0;
     Nrm=1;
     for(i=0;i<=s;i++)
      {Nrm=Nrm+(2.0*i+p)*Cm[2*i+p];}

     Nrm=1/Nrm;
    }

  if ( typ=='o' )
    {
     Nrm=Cm[p];
     for (i=1; i<=s; i++)
      { Nrm=Nrm+Cm[2*i+p];}

     Nrm=1/Nrm;
    }

    Norm=Nrm;
    delete[] w;
    return conv;

}
/*==================================================================*/
#undef h_type
SHAR_EOF
fi # end of overwriting check
if test -f 'mmathur.h'
then
	echo shar: will not over-write existing file "'mmathur.h'"
else
cat << "SHAR_EOF" > 'mmathur.h'

   double MathuMZn(char typ,int n, double h, double u, short kind);
   double MathuMZDn(char typ,int n, double h, double u, short kind);

   double MathuQn(char typ,int n, double h, double v, short kind);
   double MathuQDn(char typ,int n, double h, double v, short kind);

   double  MathuIn(char typ, int n, double Am[], double Ims[], double Zm[], int CDIM);
   double  MathuKn(char typ, int n, double Am[], double Kms[], double Zm[], int CDIM);
   double  MathuIDn(char typ, int n, double h, double u, double Am[], double Ims[], double IDms[], double Zm[], double ZDm[], int CDIM);
   double  MathuKDn(char typ, int n, double h, double u, double Am[], double Ims[], double IDms[], double Zm[], double ZDm[], int CDIM);

   double  MathuQn(char typ, int n, double v, double Am[], int CDIM);
   double  MathuQDn(char typ, int n, double v, double Am[], int CDIM);
   double  IQn(char typ,int ni, double vi, double ti, double ui, double Am[], int CDIM);
   double  MathuEn(char typ, int n, double v, double Cm[], double Norm, double Am[], int CDIM);
   double  MathuEDn(char typ, int n, double v, double Cm[], double Norm, double Am[], int CDIM);


   double MCoefficients(char typ,int n, double h, double Am[], int CDIM, double *mcnr, double *Norm);
   int MCoefficientSec(char typ, int n, double h, double mcn, double Am[], int CDIM, double Cm[], double *Norm);
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
if test ! -d 'Doc'
then
	mkdir 'Doc'
fi
cd 'Doc'
if test -f 'Read.me'
then
	echo shar: will not over-write existing file "'Read.me'"
else
cat << "SHAR_EOF" > 'Read.me'
Last update 8/5/2000
Author: Fayez A. Alhargan
email: alhargan@kacst.edu.sa
Tel: 966-1-4813770
Fax: 966-1-4813764

Programs for the papers:
Algorithms for the Computation of all Mathieu Functions of Integer Orders
and
Subroutines for the computation of  Mathieu functions of integer orders
which are published in ACM Trans. On Mathematical Software.

The programs are writen in C++ and contained in the following files:
  MATHUR.CPP     Mathieu functions
  MMATHUR.CPP    Modified Mathieu functions
  MCNR.CPP       Mathieu charactristic Numbers
  BSSLR.CPP      Bessel functions
  MBSSLR.CPP     Modified Bessel functions
Plus the header files:
  MATHUR.h     Mathieu functions
  MMATHUR.h    Modified Mathieu functions
  MCNR.h       Mathieu charactristic Numbers
  BSSLR.h      Bessel functions
  MBSSLR.h     Modified Bessel functions

The driver for the programs is
  main.cpp

It has been tested using Borland C++ v3.1 using the project file
   drive.prj

The test data input/output of the program are as follows:

  input parameters:(e,3,c,f,s)
  output file name:E03FCS00.DAT (SEE BELOW)
  contains:
  Circumferential First  kind  Even   Mathieu function of order n=3

  input parameters:(o,3,d,r,m)
  output file name:O03DRM00.DAT (SEE BELOW)
  contains:
  Modified  Radial Second  kind  Odd   Mathieu function of order n=3

NOTE:  main program changed to take data from stdin and output on stdout.
SHAR_EOF
fi # end of overwriting check
cd ..
#       End of shell archive
exit 0