Subroutine PSeries0(s,L,p,pp,p2p) c This subroutine computes the value of p using the power series c solutions for the case omega = 0. Implicit None Integer i,L,nom Real*8 p,pp,p2p,a2,a3,a4,a5,a6,a7,a8 Real*8 b2,b3,b4,b5,b6,b7,b8 Real*8 p0,p1,p2,p3,p4,p5,p6,p7 Real*8 alpha,ll1,xl,xi,s Real*8 omi,delom,test2 COMMON/PAR/xi,omi,delom,nom Common/coeff/a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8 c Note that PSeries uses a8 and b8, this program does not xl = L ll1 = xl*(xl+1.d0) test2 = b2**2 + 4*b2*ll1 + 8*b2*xi - 8*b2**2*xi If(dabs(test2).lt.1.d-12) then alpha = -0.5d0 Else alpha = (-b2 + dsqrt(b2**2 + 4*b2*ll1 + 8*b2*xi - 8*b2**2*xi))/(2.*b2) Endif c THIS PROGRAM IS SET UP FOR a2 = 1!!! CHANGE IT p0 = 1.d0 p1 = (-4*alpha*b2*p0 - a3*alpha*b2*p0 - 3*alpha*b3*p0 - - 2*alpha**2*b3*p0 - 4*ll1*p0 - 8*p0*xi - 16*b2*p0*xi - - 6*a3*b2*p0*xi - 6*b3*p0*xi)/ - (2.*(2*b2 + 3*alpha*b2 + alpha**2*b2 - ll1 - 2*xi + 2*b2*xi)) p2 = (4*alpha*b2*p0 + a3**2*alpha*b2*p0 - 2*a4*alpha*b2*p0 - - 4*alpha*b3*p0 - a3*alpha*b3*p0 - 4*alpha*b4*p0 - - 2*alpha**2*b4*p0 + 6*ll1*p0 - 4*b2*p1 - a3*b2*p1 - 4*alpha*b2*p1 - - a3*alpha*b2*p1 - 5*b3*p1 - 7*alpha*b3*p1 - 2*alpha**2*b3*p1 - - 4*ll1*p1 + 12*p0*xi + 12*b2*p0*xi - 4*a3*b2*p0*xi + - 7*a3**2*b2*p0*xi - 16*a4*b2*p0*xi - 20*b3*p0*xi - 7*a3*b3*p0*xi - -8*b4*p0*xi - 8*p1*xi - 16*b2*p1*xi - 6*a3*b2*p1*xi - 6*b3*p1*xi)/ - (2.*(6*b2 + 5*alpha*b2 + alpha**2*b2 - ll1 - 2*xi + 2*b2*xi)) p3 = (-4*alpha*b2*p0 - a3**3*alpha*b2*p0 + 3*a3*a4*alpha*b2*p0 - - 3*a5*alpha*b2*p0 + 4*alpha*b3*p0 + a3**2*alpha*b3*p0 - - 2*a4*alpha*b3*p0 - 4*alpha*b4*p0 - a3*alpha*b4*p0 - - 5*alpha*b5*p0 - 2*alpha**2*b5*p0 - 8*ll1*p0 + 4*b2*p1 + - a3**2*b2*p1 - 2*a4*b2*p1 + 4*alpha*b2*p1 + a3**2*alpha*b2*p1 - - 2*a4*alpha*b2*p1 - 4*b3*p1 - a3*b3*p1 - 4*alpha*b3*p1 - - a3*alpha*b3*p1 - 6*b4*p1 - 8*alpha*b4*p1 - 2*alpha**2*b4*p1 + - 6*ll1*p1 - 8*b2*p2 - 2*a3*b2*p2 - 4*alpha*b2*p2 - a3*alpha*b2*p2 - - 14*b3*p2 - 11*alpha*b3*p2 - 2*alpha**2*b3*p2 - 4*ll1*p2 - - 16*p0*xi - 8*b2*p0*xi + 4*a3*b2*p0*xi + 4*a3**2*b2*p0*xi - - 8*a3**3*b2*p0*xi - 8*a4*b2*p0*xi + 26*a3*a4*b2*p0*xi - - 30*a5*b2*p0*xi + 16*b3*p0*xi - 4*a3*b3*p0*xi + 8*a3**2*b3*p0*xi - -18*a4*b3*p0*xi - 24*b4*p0*xi - 8*a3*b4*p0*xi - 10*b5*p0*xi + - 12*p1*xi + 12*b2*p1*xi - 4*a3*b2*p1*xi + 7*a3**2*b2*p1*xi - - 16*a4*b2*p1*xi - 20*b3*p1*xi - 7*a3*b3*p1*xi - 8*b4*p1*xi - - 8*p2*xi - 16*b2*p2*xi - 6*a3*b2*p2*xi - 6*b3*p2*xi)/ - (2.*(12*b2 + 7*alpha*b2 + alpha**2*b2 - ll1 - 2*xi + 2*b2*xi)) p4 = (4*alpha*b2*p0 + a3**4*alpha*b2*p0 - 4*a3**2*a4*alpha*b2*p0 + - 2*a4**2*alpha*b2*p0 + 4*a3*a5*alpha*b2*p0 - 4*a6*alpha*b2*p0 - - 4*alpha*b3*p0 - a3**3*alpha*b3*p0 + 3*a3*a4*alpha*b3*p0 - - 3*a5*alpha*b3*p0 + 4*alpha*b4*p0 + a3**2*alpha*b4*p0 - - 2*a4*alpha*b4*p0 - 4*alpha*b5*p0 - a3*alpha*b5*p0 - - 6*alpha*b6*p0 - 2*alpha**2*b6*p0 + 10*ll1*p0 - 4*b2*p1 - - a3**3*b2*p1 + 3*a3*a4*b2*p1 - 3*a5*b2*p1 - 4*alpha*b2*p1 - - a3**3*alpha*b2*p1 + 3*a3*a4*alpha*b2*p1 - 3*a5*alpha*b2*p1 + - 4*b3*p1 + a3**2*b3*p1 - 2*a4*b3*p1 + 4*alpha*b3*p1 + - a3**2*alpha*b3*p1 - 2*a4*alpha*b3*p1 - 4*b4*p1 - a3*b4*p1 - - 4*alpha*b4*p1 - a3*alpha*b4*p1 - 7*b5*p1 - 9*alpha*b5*p1 - - 2*alpha**2*b5*p1 - 8*ll1*p1 + 8*b2*p2 + 2*a3**2*b2*p2 - - 4*a4*b2*p2 + 4*alpha*b2*p2 + a3**2*alpha*b2*p2 - - 2*a4*alpha*b2*p2 - 8*b3*p2 - 2*a3*b3*p2 - 4*alpha*b3*p2 - - a3*alpha*b3*p2 - 16*b4*p2 - 12*alpha*b4*p2 - 2*alpha**2*b4*p2 + - 6*ll1*p2 - 12*b2*p3 - 3*a3*b2*p3 - 4*alpha*b2*p3 - - a3*alpha*b2*p3 - 27*b3*p3 - 15*alpha*b3*p3 - 2*alpha**2*b3*p3 - - 4*ll1*p3 + 20*p0*xi + 4*b2*p0*xi - 4*a3*b2*p0*xi - - 4*a3**2*b2*p0*xi - 4*a3**3*b2*p0*xi + 9*a3**4*b2*p0*xi + - 8*a4*b2*p0*xi + 12*a3*a4*b2*p0*xi - 38*a3**2*a4*b2*p0*xi + - 20*a4**2*b2*p0*xi - 12*a5*b2*p0*xi + 42*a3*a5*b2*p0*xi - - 48*a6*b2*p0*xi - 12*b3*p0*xi + 4*a3*b3*p0*xi + 4*a3**2*b3*p0*xi - - 9*a3**3*b3*p0*xi - 8*a4*b3*p0*xi + 29*a3*a4*b3*p0*xi - - 33*a5*b3*p0*xi + 20*b4*p0*xi - 4*a3*b4*p0*xi + 9*a3**2*b4*p0*xi - -20*a4*b4*p0*xi - 28*b5*p0*xi - 9*a3*b5*p0*xi - 12*b6*p0*xi - - 16*p1*xi - 8*b2*p1*xi + 4*a3*b2*p1*xi + 4*a3**2*b2*p1*xi - - 8*a3**3*b2*p1*xi - 8*a4*b2*p1*xi + 26*a3*a4*b2*p1*xi - - 30*a5*b2*p1*xi + 16*b3*p1*xi - 4*a3*b3*p1*xi + 8*a3**2*b3*p1*xi - -18*a4*b3*p1*xi - 24*b4*p1*xi - 8*a3*b4*p1*xi - 10*b5*p1*xi + - 12*p2*xi + 12*b2*p2*xi - 4*a3*b2*p2*xi + 7*a3**2*b2*p2*xi - - 16*a4*b2*p2*xi - 20*b3*p2*xi - 7*a3*b3*p2*xi - 8*b4*p2*xi - - 8*p3*xi - 16*b2*p3*xi - 6*a3*b2*p3*xi - 6*b3*p3*xi)/ - (2.*(20*b2 + 9*alpha*b2 + alpha**2*b2 - ll1 - 2*xi + 2*b2*xi)) p5 = (-4*alpha*b2*p0 - a3**5*alpha*b2*p0 + 5*a3**3*a4*alpha*b2*p0 - - 5*a3*a4**2*alpha*b2*p0 - 5*a3**2*a5*alpha*b2*p0 + - 5*a4*a5*alpha*b2*p0 + 5*a3*a6*alpha*b2*p0 - 5*a7*alpha*b2*p0 + - 4*alpha*b3*p0 + a3**4*alpha*b3*p0 - 4*a3**2*a4*alpha*b3*p0 + - 2*a4**2*alpha*b3*p0 + 4*a3*a5*alpha*b3*p0 - 4*a6*alpha*b3*p0 - - 4*alpha*b4*p0 - a3**3*alpha*b4*p0 + 3*a3*a4*alpha*b4*p0 - - 3*a5*alpha*b4*p0 + 4*alpha*b5*p0 + a3**2*alpha*b5*p0 - - 2*a4*alpha*b5*p0 - 4*alpha*b6*p0 - a3*alpha*b6*p0 - - 7*alpha*b7*p0 - 2*alpha**2*b7*p0 - 12*ll1*p0 + 4*b2*p1 + - a3**4*b2*p1 - 4*a3**2*a4*b2*p1 + 2*a4**2*b2*p1 + 4*a3*a5*b2*p1 - - 4*a6*b2*p1 + 4*alpha*b2*p1 + a3**4*alpha*b2*p1 - - 4*a3**2*a4*alpha*b2*p1 + 2*a4**2*alpha*b2*p1 + - 4*a3*a5*alpha*b2*p1 - 4*a6*alpha*b2*p1 - 4*b3*p1 - a3**3*b3*p1 + - 3*a3*a4*b3*p1 - 3*a5*b3*p1 - 4*alpha*b3*p1 - a3**3*alpha*b3*p1 + - 3*a3*a4*alpha*b3*p1 - 3*a5*alpha*b3*p1 + 4*b4*p1 + a3**2*b4*p1 - - 2*a4*b4*p1 + 4*alpha*b4*p1 + a3**2*alpha*b4*p1 - - 2*a4*alpha*b4*p1 - 4*b5*p1 - a3*b5*p1 - 4*alpha*b5*p1 - - a3*alpha*b5*p1 - 8*b6*p1 - 10*alpha*b6*p1 - 2*alpha**2*b6*p1 + - 10*ll1*p1 - 8*b2*p2 - 2*a3**3*b2*p2 + 6*a3*a4*b2*p2 - 6*a5*b2*p2 - - 4*alpha*b2*p2 - a3**3*alpha*b2*p2 + 3*a3*a4*alpha*b2*p2 - - 3*a5*alpha*b2*p2 + 8*b3*p2 + 2*a3**2*b3*p2 - 4*a4*b3*p2 + - 4*alpha*b3*p2 + a3**2*alpha*b3*p2 - 2*a4*alpha*b3*p2 - 8*b4*p2 - - 2*a3*b4*p2 - 4*alpha*b4*p2 - a3*alpha*b4*p2 - 18*b5*p2 - - 13*alpha*b5*p2 - 2*alpha**2*b5*p2 - 8*ll1*p2 + 12*b2*p3 + - 3*a3**2*b2*p3 - 6*a4*b2*p3 + 4*alpha*b2*p3 + a3**2*alpha*b2*p3 - - 2*a4*alpha*b2*p3 - 12*b3*p3 - 3*a3*b3*p3 - 4*alpha*b3*p3 - - a3*alpha*b3*p3 - 30*b4*p3 - 16*alpha*b4*p3 - 2*alpha**2*b4*p3 + - 6*ll1*p3 - 16*b2*p4 - 4*a3*b2*p4 - 4*alpha*b2*p4 - - a3*alpha*b2*p4 - 44*b3*p4 - 19*alpha*b3*p4 - 2*alpha**2*b3*p4 - - 4*ll1*p4 - 24*p0*xi + 4*a3*b2*p0*xi + 4*a3**2*b2*p0*xi + - 4*a3**3*b2*p0*xi + 4*a3**4*b2*p0*xi - 10*a3**5*b2*p0*xi - - 8*a4*b2*p0*xi - 12*a3*a4*b2*p0*xi - 16*a3**2*a4*b2*p0*xi + - 52*a3**3*a4*b2*p0*xi + 8*a4**2*b2*p0*xi - 54*a3*a4**2*b2*p0*xi + - 12*a5*b2*p0*xi + 16*a3*a5*b2*p0*xi - 56*a3**2*a5*b2*p0*xi + - 58*a4*a5*b2*p0*xi - 16*a6*b2*p0*xi + 62*a3*a6*b2*p0*xi - - 70*a7*b2*p0*xi + 8*b3*p0*xi - 4*a3*b3*p0*xi - 4*a3**2*b3*p0*xi - - 4*a3**3*b3*p0*xi + 10*a3**4*b3*p0*xi + 8*a4*b3*p0*xi + - 12*a3*a4*b3*p0*xi - 42*a3**2*a4*b3*p0*xi + 22*a4**2*b3*p0*xi - - 12*a5*b3*p0*xi + 46*a3*a5*b3*p0*xi - 52*a6*b3*p0*xi - - 16*b4*p0*xi + 4*a3*b4*p0*xi + 4*a3**2*b4*p0*xi - - 10*a3**3*b4*p0*xi - 8*a4*b4*p0*xi + 32*a3*a4*b4*p0*xi - - 36*a5*b4*p0*xi + 24*b5*p0*xi - 4*a3*b5*p0*xi + 10*a3**2*b5*p0*xi - -22*a4*b5*p0*xi - 32*b6*p0*xi - 10*a3*b6*p0*xi - 14*b7*p0*xi + - 20*p1*xi + 4*b2*p1*xi - 4*a3*b2*p1*xi - 4*a3**2*b2*p1*xi - - 4*a3**3*b2*p1*xi + 9*a3**4*b2*p1*xi + 8*a4*b2*p1*xi + - 12*a3*a4*b2*p1*xi - 38*a3**2*a4*b2*p1*xi + 20*a4**2*b2*p1*xi - - 12*a5*b2*p1*xi + 42*a3*a5*b2*p1*xi - 48*a6*b2*p1*xi - - 12*b3*p1*xi + 4*a3*b3*p1*xi + 4*a3**2*b3*p1*xi - - 9*a3**3*b3*p1*xi - 8*a4*b3*p1*xi + 29*a3*a4*b3*p1*xi - - 33*a5*b3*p1*xi + 20*b4*p1*xi - 4*a3*b4*p1*xi + 9*a3**2*b4*p1*xi - -20*a4*b4*p1*xi - 28*b5*p1*xi - 9*a3*b5*p1*xi - 12*b6*p1*xi - - 16*p2*xi - 8*b2*p2*xi + 4*a3*b2*p2*xi + 4*a3**2*b2*p2*xi - - 8*a3**3*b2*p2*xi - 8*a4*b2*p2*xi + 26*a3*a4*b2*p2*xi - - 30*a5*b2*p2*xi + 16*b3*p2*xi - 4*a3*b3*p2*xi + 8*a3**2*b3*p2*xi - -18*a4*b3*p2*xi - 24*b4*p2*xi - 8*a3*b4*p2*xi - 10*b5*p2*xi + - 12*p3*xi + 12*b2*p3*xi - 4*a3*b2*p3*xi + 7*a3**2*b2*p3*xi - - 16*a4*b2*p3*xi - 20*b3*p3*xi - 7*a3*b3*p3*xi - 8*b4*p3*xi - - 8*p4*xi - 16*b2*p4*xi - 6*a3*b2*p4*xi - 6*b3*p4*xi)/ - (2.*(30*b2 + 11*alpha*b2 + alpha**2*b2 - ll1 - 2*xi + 2*b2*xi)) p = s**alpha*(p0+p1*s+p2*s**2+p3*s**3+p4*s**4+p5*s**5) pp = alpha*s**(alpha-1.d0)*(p0) pp = alpha*s**(alpha-1.d0)*(p0+p1*s+p2*s**2+p3*s**3+p4*s**4+p5*s**5) - + s**alpha*(p1+2.*p2*s+3.*p3*s**2+4.*p4*s**3+5.*p5*s**4) p2p = alpha*(alpha-1.d0)*s**(alpha-2.d0)*(p0+p1*s+p2*s**2+p3*s**3+p4*s**4+p5*s**5) - + 2.*alpha*s**(alpha-1.d0)*(p1+2.*p2*s+3.*p3*s**2+4.*p4*s**3+5.*p5*s**4) - + s**alpha*(2.*p2+6.*p3*s+12.*p4*s**2+20.*p5*s**3) c Print *,'as are: ',a2,a3,a4,a5,a6,a7,a8 c Print *,'bs are: ',b2,b3,b4,b5,b6,b7,b8 c Print *,'alpha = ',alpha,' p0 = ',p0,' p1 = ',p1,' p2 = ',p2 c Print*,'p3 = ',p3,' p4 = ',p4,' p5 = ',p5 c Print*,'p = ',p,' pp = ',pp,' p2p = ',p2p Return End