c
c
      program dOPairy
c
c
      external dquad
      double precision x0,x1,a1(160),b1(160),a2(160),
     *b2(160),e(160),d1mach,depsma,dendl(3),dendr(3),
     *deps,dxp(1),dyp(1),dxfer(160),dwfer(160),
     *dalpha(40),dbeta(40),dbe(40),dx(160),dw(160),
     *dxm(480),dwm(480),dp0(480),dp1(480),dp2(480)
      logical finld,finrd
      common/d/x0,x1,a1,b1,a2,b2,e
c
c This routine generates in double precision the first 40
c recurrence coefficients of the orthogonal polynomials for 
c the Airy weight function. Printed are the double-precision
c values of the alpha- and beta-coefficients. The routine
c requires the ORTHPOL routines d1mach, drecur, dgauss, dmcdis,
c dqgp, dsti, and dlancz (see W. Gautschi, "Algorithm 726: 
c ORTHPOL -- a package of routines for generating orthogonal
c polynomials and Gauss-type quadrature rules", ACM Trans. Math. 
c Software, v. 20 (1994), pp. 21-62).
c
      nmax=40
      x0=2.d0
      x1=6.d0
      call drecur(160,1,0.d0,0.d0,a1,b1,ier1)
      open(unit=5,file='dablee',status='old')
      read(*,1) (a2(k),b2(k),k=1,160)
    1 format(1x,2d23.15)
      finld=.true.
      finrd=.false.
      depsma=d1mach(3)
c
c depsma is the machine double precision.
c
      iq=1
      idelta=1
      irout=1
      mcd=3
      mp=0
      ncpmd=160
c
c Set up the partition for the discretization of the inner product.
c
      dendl(1)=0.d0
      dendr(1)=2.d0
      dendl(2)=2.d0
      dendr(2)=6.d0
      dendl(3)=6.d0
      deps=1.d3*depsma
c
c Compute the desired recursion coefficients by the multiple-component
c discretization procedure. 
c
      call dmcdis(nmax,ncpmd,mcd,mp,dxp,dyp,dquad,deps,iq,idelta,irout,
     *  finld,finrd,dendl,dendr,dxfer,dwfer,dalpha,dbeta,ncapd,kountd,
     *  ierrd,ied,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2)
c
c Print the results.
c
      write(*,2)
    2 format(/5x,'k',7x,'dalpha(k)',14x,'dbeta(k)'/)
      do 20 k=1,nmax
        km1=k-1
        write(*,3) km1,dalpha(k),dbeta(k)
    3   format(1x,i5,2d22.14)
   20 continue
      stop
      end 

      subroutine dquad(n,dx,dw,i,ierr)
      double precision dx(n),dw(n),x0,x1,a1(160),b1(160),
     *a2(160),b2(160),e(160),z(22),w(22),pi,c,c1,dgamma,
     *eps,d1mach,x,dairy
      common/d/x0,x1,a1,b1,a2,b2,e
      common/da/z,w
      pi=4.d0*datan(1.d0)
      c=.5d0**(2.d0/3.d0)
      c1=(2.d0**(2.d0/3.d0))*pi/((3.d0**(1.d0/6.d0))*dgamma(
     *  2.d0/3.d0,ie))
      eps=d1mach(3)
      if(i.eq.1) then
        call dgauss(n,a1,b1,eps,dx,dw,ierg,e)
        do 10 k=1,n
          x=.5d0*x0*(1.d0+dx(k))
          dw(k)=.5d0*x0*c1*dairy(c*(x**2))*dexp(
     *      -(x**3)/3.d0)*dw(k)
          dx(k)=x**3/3.d0
   10   continue
      else if(i.eq.2) then
        call dgauss(n,a1,b1,eps,dx,dw,ierg,e)
        do 20 k=1,n
          x=.5d0*((x1-x0)*dx(k)+x1+x0)
          dw(k)=.5d0*(x1-x0)*c1*dairy(c*(x**2))*dexp(
     *      -(x**3)/3.d0)*dw(k)
          dx(k)=x**3/3.d0
   20   continue
      else
        call dgauss(n,a2,b2,eps,dx,dw,ierg,e)
        do 30 k=1,n
          x=x1+dx(k)
          dw(k)=c1*dexp(-(x1**3)/3.d0)*dairy(c*(x**2))*
     *      dexp(-x1*(x1*dx(k)+dx(k)**2))*dw(k)
          dx(k)=(x**3)/3.d0
   30 continue
      end if
      return
      end

      double precision function dairy(x)
      double precision z(22),w(22),x,x1,pi,c1,c2,c3,
     *series,dgamma,sum
      common/da/z,w
      data z/0.511887709830339D-01,0.308510111483174D+00,
     *0.787420768335179D+00,0.149044228440672D+01,
     *0.242118226396883D+01,0.358453679177000D+01,
     *0.498684812051623D+01,0.663611185696245D+01,
     *0.854225679370883D+01,0.107175242658754D+02,
     *0.131769867508666D+02,0.159392671752419D+02,
     *0.190275572217344D+02,0.224710978067566D+02,
     *0.263074050499585D+02,0.305857613946155D+02,
     *0.353729920062386D+02,0.407637095509569D+02,
     *0.469002619672993D+02,0.540170851713205D+02,
     *0.625619516130729D+02,0.736832355972619D+02/
      data w/0.229628584688356D+00,0.328630853320906D+00,
     *0.279589369270967D+00,0.172031321598539D+00,
     *0.801806232510865D-01,0.287486894867489D-01,
     *0.796171780426947D-02,0.169980517736777D-02,
     *0.278088482408926D-03,0.345325854146748D-04,
     *0.321279872294375D-05,0.220181421137924D-06,
     *0.108763443450499D-07,0.376618119364524D-09,
     *0.881839196298758D-11,0.133138663606891D-12,
     *0.121510124790229D-14,0.611861104951773D-17,
     *0.148223889569632D-19,0.138011407375144D-22,
     *0.322948437472136D-26,0.651796266021172D-31/
      if(x.le.0.d0) then
        write(*,1)
    1   format(1x,'invalid x-argument in dairy')
        stop
      end if
      pi=4.d0*datan(1.d0)
      c1=1.d0/((3.d0**(2.d0/3.d0))*dgamma(2.d0/3.d0,ie))
      c2=1.d0/((3.d0**(4.d0/3.d0))*dgamma(4.d0/3.d0,ie))
      if(x.le.3.d0**(2.d0/3.d0)) then
        dairy=c1*series(x,-1.d0/3.d0)-c2*x*series(x,1.d0/3.d0)
        return
      else
        n=22
        c3=1.d0/(dsqrt(pi)*(48.d0**(1.d0/6.d0))*dgamma(5.d0/6.d0,ie))
        sum=0.d0
        x1=2.d0*((dsqrt(x))**3)/3.d0
        do 10 k=1,n
          sum=sum+w(k)*(2.d0+z(k)/x1)**(-1.d0/6.d0)
   10   continue
        dairy=c3*dexp(-x1)*(x1**(-1.d0/6.d0))*sum
      end if
      return
      end

      double precision function series(x,dnu)
      double precision x,dnu,x3,sum1,t1,dk,sum,t
      k=0
      x3=x**3
      sum1=1.d0
      t1=1.d0
   10 k=k+1
      dk=dble(k)
      sum=sum1
      t=t1
      t1=x3*t/(9.d0*dk*(dnu+dk))
      sum1=sum+t1
      if(sum1.ne.sum) goto 10
      series=sum1
      return
      end
  
      double precision function dwf(x,i)
      double precision x
      return
      end
c
c
                      file dablee 
         (containing the recursion coefficients
       for the weight function  w(x)=exp(-x**3/3)
c
c
      0.729011132947227D+00  0.128789931685407D+01
      0.104221982567474D+01  0.245000979417421D+00
      0.125373064220196D+01  0.353073517279907D+00
      0.140618208893400D+01  0.453806544754720D+00
      0.153047170886983D+01  0.546709151632936D+00
      0.163711468769310D+01  0.632791431265656D+00
      0.173132652800093D+01  0.713591550241559D+00
      0.181621572840940D+01  0.790171600818179D+00
      0.189380331629451D+01  0.863276695500399D+00
      0.196548682633124D+01  0.933452973983708D+00
      0.203227833945824D+01  0.100111430192640D+01
      0.209493741056696D+01  0.106658307630641D+01
      0.215405051280269D+01  0.113011630551701D+01
      0.221008111612034D+01  0.119192288798298D+01
      0.226340263870694D+01  0.125217543914883D+01
      0.231432091013042D+01  0.131101862259663D+01
      0.236308994699168D+01  0.136857516761462D+01
      0.240992331675625D+01  0.142495031218022D+01
      0.245500250496030D+01  0.148023514736833D+01
      0.249848319498193D+01  0.153450917822770D+01
      0.254050006122831D+01  0.158784231475101D+01
      0.258117048278474D+01  0.164029644098122D+01
      0.262059745939295D+01  0.169192666687487D+01
      0.265887192879668D+01  0.174278233817534D+01
      0.269607462849073D+01  0.179290785929909D+01
      0.273227760630482D+01  0.184234337002653D+01
      0.276754545717175D+01  0.189112530665539D+01
      0.280193634412483D+01  0.193928687094024D+01
      0.283550284760913D+01  0.198685842476269D+01
      0.286829267695869D+01  0.203386782448106D+01
      0.290034927029989D+01  0.208034070590724D+01
      0.293171230344481D+01  0.212630072857875D+01
      0.296241812401779D+01  0.217176978624682D+01
      0.299250012375155D+01  0.221676818914791D+01
      0.302198905933231D+01  0.226131482257079D+01
      0.305091333018266D+01  0.230542728539943D+01
      0.307929922000566D+01  0.234912201165387D+01
      0.310717110767653D+01  0.239241437752433D+01
      0.313455165208123D+01  0.243531879597110D+01
      0.316146195471043D+01  0.247784880062069D+01
      0.318792170317849D+01  0.252001712040978D+01
      0.321394929831879D+01  0.256183574620110D+01
      0.323956196708347D+01  0.260331599040762D+01
      0.326477586312814D+01  0.264446854050621D+01
      0.328960615667592D+01  0.268530350719359D+01
      0.331406711501741D+01  0.272583046782938D+01
      0.333817217480591D+01  0.276605850572120D+01
      0.336193400714172D+01  0.280599624573085D+01
      0.338536457630075D+01  0.284565188661641D+01
      0.340847519284574D+01  0.288503323047080D+01
      0.343127656175960D+01  0.292414770957096D+01
      0.345377882615615D+01  0.296300241091235D+01
      0.347599160705250D+01  0.300160409866944D+01
      0.349792403962555D+01  0.303995923479377D+01
      0.351958480632325D+01  0.307807399793612D+01
      0.354098216715593D+01  0.311595430085732D+01
      0.356212398745396D+01  0.315360580647365D+01
      0.358301776334475D+01  0.319103394266626D+01
      0.360367064517267D+01  0.322824391596956D+01
      0.362408945906008D+01  0.326524072424137D+01
      0.364428072678576D+01  0.330202916840641D+01
      0.366425068413760D+01  0.333861386335501D+01
      0.368400529787925D+01  0.337499924807086D+01
      0.370355028145616D+01  0.341118959505352D+01
      0.372289110955262D+01  0.344718901909540D+01
      0.374203303160062D+01  0.348300148546671D+01
      0.376098108433051D+01  0.351863081755667D+01
      0.377974010344498D+01  0.355408070401491D+01
      0.379831473448954D+01  0.358935470543261D+01
      0.381670944298562D+01  0.362445626059935D+01
      0.383492852388629D+01  0.365938869236831D+01
      0.385297611040875D+01  0.369415521315971D+01
      0.387085618229271D+01  0.372875893012936D+01
      0.388857257352953D+01  0.376320285002715D+01
      0.390612897960254D+01  0.379748988376813D+01
      0.392352896427572D+01  0.383162285073666D+01
      0.394077596596450D+01  0.386560448284268D+01
      0.395787330371940D+01  0.389943742834739D+01
      0.397482418285084D+01  0.393312425547423D+01
      0.399163170022085D+01  0.396666745581992D+01
      0.400829884922534D+01  0.400006944757885D+01
      0.402482852448863D+01  0.403333257859331D+01
      0.404122352629026D+01  0.406645912924110D+01
      0.405748656474219D+01  0.409945131517076D+01
      0.407362026373355D+01  0.413231128989461D+01
      0.408962716465820D+01  0.416504114724819D+01
      0.410550972993966D+01  0.419764292372475D+01
      0.412127034636651D+01  0.423011860069237D+01
      0.413691132825050D+01  0.426247010650092D+01
      0.415243492041879D+01  0.429469931848560D+01
      0.416784330105062D+01  0.432680806487306D+01
      0.418313858436826D+01  0.435879812659606D+01
      0.419832282319118D+01  0.439067123902183D+01
      0.421339801136175D+01  0.442242909359937D+01
      0.422836608605034D+01  0.445407333942998D+01
      0.424322892994689D+01  0.448560558476568D+01
      0.425798837334578D+01  0.451702739843944D+01
      0.427264619613020D+01  0.454834031123091D+01
      0.428720412966183D+01  0.457954581717136D+01
      0.430166385858128D+01  0.461064537479101D+01
      0.431602702252432D+01  0.464164040831186D+01
      0.433029521775877D+01  0.467253230878897D+01
      0.434446999874622D+01  0.470332243520288D+01
      0.435855287963296D+01  0.473401211550570D+01
      0.437254533567386D+01  0.476460264762327D+01
      0.438644880459286D+01  0.479509530041570D+01
      0.440026468788342D+01  0.482549131459827D+01
      0.441399435205224D+01  0.485579190362485D+01
      0.442763912980905D+01  0.488599825453550D+01
      0.444120032120544D+01  0.491611152877025D+01
      0.445467919472526D+01  0.494613286295051D+01
      0.446807698832912D+01  0.497606336962976D+01
      0.448139491045530D+01  0.500590413801503D+01
      0.449463414097919D+01  0.503565623466048D+01
      0.450779583213353D+01  0.506532070413446D+01
      0.452088110939111D+01  0.509489856966120D+01
      0.453389107231195D+01  0.512439083373848D+01
      0.454682679535668D+01  0.515379847873214D+01
      0.455968932866763D+01  0.518312246744866D+01
      0.457247969881929D+01  0.521236374368682D+01
      0.458519890953954D+01  0.524152323276917D+01
      0.459784794240297D+01  0.527060184205444D+01
      0.461042775749771D+01  0.529960046143163D+01
      0.462293929406687D+01  0.532851996379651D+01
      0.463538347112584D+01  0.535736120551141D+01
      0.464776118805650D+01  0.538612502684899D+01
      0.466007332517933D+01  0.541481225242060D+01
      0.467232074430460D+01  0.544342369158996D+01
      0.468450428926331D+01  0.547196013887274D+01
      0.469662478641892D+01  0.550042237432264D+01
      0.470868304516078D+01  0.552881116390451D+01
      0.472067985837982D+01  0.555712725985508D+01
      0.473261600292749D+01  0.558537140103173D+01
      0.474449224005854D+01  0.561354431324986D+01
      0.475630931585839D+01  0.564164670960928D+01
      0.476806796165564D+01  0.566967929081008D+01
      0.477976889442050D+01  0.569764274545831D+01
      0.479141281714956D+01  0.572553775036202D+01
      0.480300041923754D+01  0.575336497081790D+01
      0.481453237683656D+01  0.578112506088893D+01
      0.482600935320344D+01  0.580881866367339D+01
      0.483743199903543D+01  0.583644641156560D+01
      0.484880095279495D+01  0.586400892650859D+01
      0.486011684102368D+01  0.589150682023911D+01
      0.487138027864647D+01  0.591894069452523D+01
      0.488259186926550D+01  0.594631114139677D+01
      0.489375220544494D+01  0.597361874336888D+01
      0.490486186898663D+01  0.600086407365899D+01
      0.491592143119704D+01  0.602804769639735D+01
      0.492693145314585D+01  0.605517016683144D+01
      0.493789248591645D+01  0.608223203152444D+01
      0.494880507084879D+01  0.610923382854797D+01
      0.495966973977463D+01  0.613617608766927D+01
      0.497048701524569D+01  0.616305933053315D+01
      0.498125741075486D+01  0.618988407083865D+01
      0.499198143095076D+01  0.621665081451088D+01
      0.500265957184580D+01  0.624336005986795D+01
      0.501329232101818D+01  0.627001229778335D+01
      0.502388015780777D+01  0.629660801184378D+01
      0.503442355350636D+01  0.632314767850273D+01


