Subroutine PSeries(s,omega,L,p,p1p,p2p) c This subroutine computes the value of p using the power series c solutions for the case omega not equal to zero. Implicit None Integer i,L,nom Real*8 om,p,p1p,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,omega COMMON/PAR/xi,omi,delom,nom Common/coeff/a2,a3,a4,a5,a6,a7,a8,b2,b3,b4,b5,b6,b7,b8 c c THIS PROGRAM IS SET UP FOR a2 = 1. CHANGE IT!!! c c Note that in the mathematica program that generated these series c I scaled omega so that om_new = r0^2*om_old/Sqrt[b2] c However, I don't want to do that scaling in the main program I decided c so I will scale back the b2 part. om = omega/dsqrt(b2) xl = L ll1 = xl*(xl+1.d0) alpha = -((a3*b2 + b3)*om)/(2.*b2) p0 = 1.d0 p1 = (4*b2*ll1*p0 - 8*b2**2*om*p0 + 3*a3**2*b2**2*om**2*p0 - - 4*a4*b2**2*om**2*p0 + 2*a3*b2*b3*om**2*p0 + 3*b3**2*om**2*p0 - - 4*b2*b4*om**2*p0 + 8*b2*p0*xi - 8*b2**2*p0*xi)/(8.*b2**2*om) p2 = (-8*b2**2*ll1*p0 + 8*b2**3*om*p0 + 4*a3*b2**3*om*p0 + - 3*a3**2*b2**3*om*p0 - 4*a4*b2**3*om*p0 - 4*b2**2*b3*om*p0 + - 2*a3*b2**2*b3*om*p0 + 3*b2*b3**2*om*p0 - 4*b2**2*b4*om*p0 - - 4*a3**3*b2**3*om**2*p0 + 8*a3*a4*b2**3*om**2*p0 - - 4*a5*b2**3*om**2*p0 - a3**2*b2**2*b3*om**2*p0 - - 2*a3*b2*b3**2*om**2*p0 - b3**3*om**2*p0 + 4*a3*b2**2*b4*om**2*p0 + - 4*b2*b3*b4*om**2*p0 - 4*b2**2*b5*om**2*p0 - 8*b2**3*p1 + - 4*b2**2*ll1*p1 - 8*b2**3*om*p1 + 4*a3*b2**3*om*p1 - - 4*b2**2*b3*om*p1 + 3*a3**2*b2**3*om**2*p1 - 4*a4*b2**3*om**2*p1 + - 2*a3*b2**2*b3*om**2*p1 + 3*b2*b3**2*om**2*p1 - - 4*b2**2*b4*om**2*p1 - 16*b2**2*p0*xi - 32*b2**3*p0*xi - - 12*a3*b2**3*p0*xi - 12*b2**2*b3*p0*xi + 8*b2**2*p1*xi - - 8*b2**3*p1*xi)/(16.*b2**3*om) p3 = (12*b2**2*ll1*p0 - 8*b2**3*om*p0 - 4*a3*b2**3*om*p0 - - 3*a3**3*b2**3*om*p0 + 8*a3*a4*b2**3*om*p0 - 6*a5*b2**3*om*p0 + - 4*b2**2*b3*om*p0 + 4*a3*b2**2*b3*om*p0 + 2*a3**2*b2**2*b3*om*p0 - - 2*a4*b2**2*b3*om*p0 + 4*b2*b3**2*om*p0 + a3*b2*b3**2*om*p0 - - 8*b2**2*b4*om*p0 + 2*a3*b2**2*b4*om*p0 + 4*b2*b3*b4*om*p0 - - 6*b2**2*b5*om*p0 + 4*a3**4*b2**3*om**2*p0 - - 12*a3**2*a4*b2**3*om**2*p0 + 4*a4**2*b2**3*om**2*p0 + - 8*a3*a5*b2**3*om**2*p0 - 4*a6*b2**3*om**2*p0 - - a3**2*b2**2*b4*om**2*p0 - 2*a3*b2*b3*b4*om**2*p0 - - b3**2*b4*om**2*p0 + 4*a3*b2**2*b5*om**2*p0 + 4*b2*b3*b5*om**2*p0 - - 4*b2**2*b6*om**2*p0 - 8*b2**3*p1 - 2*a3*b2**3*p1 - - 10*b2**2*b3*p1 - 8*b2**2*ll1*p1 + 8*b2**3*om*p1 + - 4*a3*b2**3*om*p1 + 3*a3**2*b2**3*om*p1 - 4*a4*b2**3*om*p1 - - 4*b2**2*b3*om*p1 + 6*a3*b2**2*b3*om*p1 + 7*b2*b3**2*om*p1 - - 12*b2**2*b4*om*p1 - 4*a3**3*b2**3*om**2*p1 + - 8*a3*a4*b2**3*om**2*p1 - 4*a5*b2**3*om**2*p1 - - a3**2*b2**2*b3*om**2*p1 - 2*a3*b2*b3**2*om**2*p1 - - b3**3*om**2*p1 + 4*a3*b2**2*b4*om**2*p1 + 4*b2*b3*b4*om**2*p1 - - 4*b2**2*b5*om**2*p1 - 24*b2**3*p2 + 4*b2**2*ll1*p2 - - 8*b2**3*om*p2 + 8*a3*b2**3*om*p2 - 8*b2**2*b3*om*p2 + - 3*a3**2*b2**3*om**2*p2 - 4*a4*b2**3*om**2*p2 + - 2*a3*b2**2*b3*om**2*p2 + 3*b2*b3**2*om**2*p2 - - 4*b2**2*b4*om**2*p2 + 24*b2**2*p0*xi + 24*b2**3*p0*xi - - 8*a3*b2**3*p0*xi + 14*a3**2*b2**3*p0*xi - 32*a4*b2**3*p0*xi - - 40*b2**2*b3*p0*xi - 14*a3*b2**2*b3*p0*xi - 16*b2**2*b4*p0*xi - - 16*b2**2*p1*xi - 32*b2**3*p1*xi - 12*a3*b2**3*p1*xi - - 12*b2**2*b3*p1*xi + 8*b2**2*p2*xi - 8*b2**3*p2*xi)/(24.*b2**3*om) p4 = (-16*b2**2*ll1*p0 + 8*b2**3*om*p0 + 4*a3*b2**3*om*p0 + - 3*a3**4*b2**3*om*p0 - 11*a3**2*a4*b2**3*om*p0 + - 4*a4**2*b2**3*om*p0 + 11*a3*a5*b2**3*om*p0 - 8*a6*b2**3*om*p0 - - 4*b2**2*b3*om*p0 - 4*a3*b2**2*b3*om*p0 - 2*a3**3*b2**2*b3*om*p0 + - 5*a3*a4*b2**2*b3*om*p0 - 3*a5*b2**2*b3*om*p0 - 4*b2*b3**2*om*p0 - - a3**2*b2*b3**2*om*p0 + 2*a4*b2*b3**2*om*p0 + 8*b2**2*b4*om*p0 + - 4*a3*b2**2*b4*om*p0 + 3*a3**2*b2**2*b4*om*p0 - - 4*a4*b2**2*b4*om*p0 + 4*b2*b3*b4*om*p0 + a3*b2*b3*b4*om*p0 - - 8*b2**2*b5*om*p0 + 3*a3*b2**2*b5*om*p0 + 5*b2*b3*b5*om*p0 - - 8*b2**2*b6*om*p0 - 4*a3**5*b2**3*om**2*p0 + - 16*a3**3*a4*b2**3*om**2*p0 - 12*a3*a4**2*b2**3*om**2*p0 - - 12*a3**2*a5*b2**3*om**2*p0 + 8*a4*a5*b2**3*om**2*p0 + - 8*a3*a6*b2**3*om**2*p0 - 4*a7*b2**3*om**2*p0 - - a3**2*b2**2*b5*om**2*p0 - 2*a3*b2*b3*b5*om**2*p0 - - b3**2*b5*om**2*p0 + 4*a3*b2**2*b6*om**2*p0 + 4*b2*b3*b6*om**2*p0 - - 4*b2**2*b7*om**2*p0 + 8*b2**3*p1 + 2*a3**2*b2**3*p1 - - 4*a4*b2**3*p1 - 8*b2**2*b3*p1 - 2*a3*b2**2*b3*p1 - - 12*b2**2*b4*p1 + 12*b2**2*ll1*p1 - 8*b2**3*om*p1 - - 4*a3*b2**3*om*p1 - 3*a3**3*b2**3*om*p1 + 8*a3*a4*b2**3*om*p1 - - 6*a5*b2**3*om*p1 + 4*b2**2*b3*om*p1 + 4*a3*b2**2*b3*om*p1 + - 2*a3**2*b2**2*b3*om*p1 - 2*a4*b2**2*b3*om*p1 + 4*b2*b3**2*om*p1 + - a3*b2*b3**2*om*p1 - 8*b2**2*b4*om*p1 + 6*a3*b2**2*b4*om*p1 + - 8*b2*b3*b4*om*p1 - 14*b2**2*b5*om*p1 + 4*a3**4*b2**3*om**2*p1 - - 12*a3**2*a4*b2**3*om**2*p1 + 4*a4**2*b2**3*om**2*p1 + - 8*a3*a5*b2**3*om**2*p1 - 4*a6*b2**3*om**2*p1 - - a3**2*b2**2*b4*om**2*p1 - 2*a3*b2*b3*b4*om**2*p1 - - b3**2*b4*om**2*p1 + 4*a3*b2**2*b5*om**2*p1 + 4*b2*b3*b5*om**2*p1 - - 4*b2**2*b6*om**2*p1 - 16*b2**3*p2 - 4*a3*b2**3*p2 - - 28*b2**2*b3*p2 - 8*b2**2*ll1*p2 + 8*b2**3*om*p2 + - 4*a3*b2**3*om*p2 + 3*a3**2*b2**3*om*p2 - 4*a4*b2**3*om*p2 - - 4*b2**2*b3*om*p2 + 10*a3*b2**2*b3*om*p2 + 11*b2*b3**2*om*p2 - - 20*b2**2*b4*om*p2 - 4*a3**3*b2**3*om**2*p2 + - 8*a3*a4*b2**3*om**2*p2 - 4*a5*b2**3*om**2*p2 - - a3**2*b2**2*b3*om**2*p2 - 2*a3*b2*b3**2*om**2*p2 - - b3**3*om**2*p2 + 4*a3*b2**2*b4*om**2*p2 + 4*b2*b3*b4*om**2*p2 - - 4*b2**2*b5*om**2*p2 - 48*b2**3*p3 + 4*b2**2*ll1*p3 - - 8*b2**3*om*p3 + 12*a3*b2**3*om*p3 - 12*b2**2*b3*om*p3 + - 3*a3**2*b2**3*om**2*p3 - 4*a4*b2**3*om**2*p3 + - 2*a3*b2**2*b3*om**2*p3 + 3*b2*b3**2*om**2*p3 - - 4*b2**2*b4*om**2*p3 - 32*b2**2*p0*xi - 16*b2**3*p0*xi + - 8*a3*b2**3*p0*xi + 8*a3**2*b2**3*p0*xi - 16*a3**3*b2**3*p0*xi - - 16*a4*b2**3*p0*xi + 52*a3*a4*b2**3*p0*xi - 60*a5*b2**3*p0*xi + - 32*b2**2*b3*p0*xi - 8*a3*b2**2*b3*p0*xi + - 16*a3**2*b2**2*b3*p0*xi - 36*a4*b2**2*b3*p0*xi - - 48*b2**2*b4*p0*xi - 16*a3*b2**2*b4*p0*xi - 20*b2**2*b5*p0*xi + - 24*b2**2*p1*xi + 24*b2**3*p1*xi - 8*a3*b2**3*p1*xi + - 14*a3**2*b2**3*p1*xi - 32*a4*b2**3*p1*xi - 40*b2**2*b3*p1*xi - - 14*a3*b2**2*b3*p1*xi - 16*b2**2*b4*p1*xi - 16*b2**2*p2*xi - - 32*b2**3*p2*xi - 12*a3*b2**3*p2*xi - 12*b2**2*b3*p2*xi + - 8*b2**2*p3*xi - 8*b2**3*p3*xi)/(32.*b2**3*om) p5 = (20*b2**2*ll1*p0 - 8*b2**3*om*p0 - 4*a3*b2**3*om*p0 - - 3*a3**5*b2**3*om*p0 + 14*a3**3*a4*b2**3*om*p0 - - 12*a3*a4**2*b2**3*om*p0 - 14*a3**2*a5*b2**3*om*p0 + - 10*a4*a5*b2**3*om*p0 + 14*a3*a6*b2**3*om*p0 - 10*a7*b2**3*om*p0 + - 4*b2**2*b3*om*p0 + 4*a3*b2**2*b3*om*p0 + 2*a3**4*b2**2*b3*om*p0 - - 7*a3**2*a4*b2**2*b3*om*p0 + 2*a4**2*b2**2*b3*om*p0 + - 7*a3*a5*b2**2*b3*om*p0 - 4*a6*b2**2*b3*om*p0 + 4*b2*b3**2*om*p0 + - a3**3*b2*b3**2*om*p0 - 3*a3*a4*b2*b3**2*om*p0 + - 3*a5*b2*b3**2*om*p0 - 8*b2**2*b4*om*p0 - 4*a3*b2**2*b4*om*p0 - - 3*a3**3*b2**2*b4*om*p0 + 8*a3*a4*b2**2*b4*om*p0 - - 6*a5*b2**2*b4*om*p0 - 4*b2*b3*b4*om*p0 - a3**2*b2*b3*b4*om*p0 + - 2*a4*b2*b3*b4*om*p0 + 8*b2**2*b5*om*p0 + 4*a3*b2**2*b5*om*p0 + - 3*a3**2*b2**2*b5*om*p0 - 4*a4*b2**2*b5*om*p0 + 4*b2*b3*b5*om*p0 + - a3*b2*b3*b5*om*p0 - 8*b2**2*b6*om*p0 + 4*a3*b2**2*b6*om*p0 + - 6*b2*b3*b6*om*p0 - 10*b2**2*b7*om*p0 + 4*a3**6*b2**3*om**2*p0 - - 20*a3**4*a4*b2**3*om**2*p0 + 24*a3**2*a4**2*b2**3*om**2*p0 - - 4*a4**3*b2**3*om**2*p0 + 16*a3**3*a5*b2**3*om**2*p0 - - 24*a3*a4*a5*b2**3*om**2*p0 + 4*a5**2*b2**3*om**2*p0 - - 12*a3**2*a6*b2**3*om**2*p0 + 8*a4*a6*b2**3*om**2*p0 + - 8*a3*a7*b2**3*om**2*p0 - 4*a8*b2**3*om**2*p0 - - a3**2*b2**2*b6*om**2*p0 - 2*a3*b2*b3*b6*om**2*p0 - - b3**2*b6*om**2*p0 + 4*a3*b2**2*b7*om**2*p0 + 4*b2*b3*b7*om**2*p0 - - 4*b2**2*b8*om**2*p0 - 8*b2**3*p1 - 2*a3**3*b2**3*p1 + - 6*a3*a4*b2**3*p1 - 6*a5*b2**3*p1 + 8*b2**2*b3*p1 + - 2*a3**2*b2**2*b3*p1 - 4*a4*b2**2*b3*p1 - 8*b2**2*b4*p1 - - 2*a3*b2**2*b4*p1 - 14*b2**2*b5*p1 - 16*b2**2*ll1*p1 + - 8*b2**3*om*p1 + 4*a3*b2**3*om*p1 + 3*a3**4*b2**3*om*p1 - - 11*a3**2*a4*b2**3*om*p1 + 4*a4**2*b2**3*om*p1 + - 11*a3*a5*b2**3*om*p1 - 8*a6*b2**3*om*p1 - 4*b2**2*b3*om*p1 - - 4*a3*b2**2*b3*om*p1 - 2*a3**3*b2**2*b3*om*p1 + - 5*a3*a4*b2**2*b3*om*p1 - 3*a5*b2**2*b3*om*p1 - 4*b2*b3**2*om*p1 - - a3**2*b2*b3**2*om*p1 + 2*a4*b2*b3**2*om*p1 + 8*b2**2*b4*om*p1 + - 4*a3*b2**2*b4*om*p1 + 3*a3**2*b2**2*b4*om*p1 - - 4*a4*b2**2*b4*om*p1 + 4*b2*b3*b4*om*p1 + a3*b2*b3*b4*om*p1 - - 8*b2**2*b5*om*p1 + 7*a3*b2**2*b5*om*p1 + 9*b2*b3*b5*om*p1 - - 16*b2**2*b6*om*p1 - 4*a3**5*b2**3*om**2*p1 + - 16*a3**3*a4*b2**3*om**2*p1 - 12*a3*a4**2*b2**3*om**2*p1 - - 12*a3**2*a5*b2**3*om**2*p1 + 8*a4*a5*b2**3*om**2*p1 + - 8*a3*a6*b2**3*om**2*p1 - 4*a7*b2**3*om**2*p1 - - a3**2*b2**2*b5*om**2*p1 - 2*a3*b2*b3*b5*om**2*p1 - - b3**2*b5*om**2*p1 + 4*a3*b2**2*b6*om**2*p1 + 4*b2*b3*b6*om**2*p1 - - 4*b2**2*b7*om**2*p1 + 16*b2**3*p2 + 4*a3**2*b2**3*p2 - - 8*a4*b2**3*p2 - 16*b2**2*b3*p2 - 4*a3*b2**2*b3*p2 - - 32*b2**2*b4*p2 + 12*b2**2*ll1*p2 - 8*b2**3*om*p2 - - 4*a3*b2**3*om*p2 - 3*a3**3*b2**3*om*p2 + 8*a3*a4*b2**3*om*p2 - - 6*a5*b2**3*om*p2 + 4*b2**2*b3*om*p2 + 4*a3*b2**2*b3*om*p2 + - 2*a3**2*b2**2*b3*om*p2 - 2*a4*b2**2*b3*om*p2 + 4*b2*b3**2*om*p2 + - a3*b2*b3**2*om*p2 - 8*b2**2*b4*om*p2 + 10*a3*b2**2*b4*om*p2 + - 12*b2*b3*b4*om*p2 - 22*b2**2*b5*om*p2 + 4*a3**4*b2**3*om**2*p2 - - 12*a3**2*a4*b2**3*om**2*p2 + 4*a4**2*b2**3*om**2*p2 + - 8*a3*a5*b2**3*om**2*p2 - 4*a6*b2**3*om**2*p2 - - a3**2*b2**2*b4*om**2*p2 - 2*a3*b2*b3*b4*om**2*p2 - - b3**2*b4*om**2*p2 + 4*a3*b2**2*b5*om**2*p2 + 4*b2*b3*b5*om**2*p2 - - 4*b2**2*b6*om**2*p2 - 24*b2**3*p3 - 6*a3*b2**3*p3 - - 54*b2**2*b3*p3 - 8*b2**2*ll1*p3 + 8*b2**3*om*p3 + - 4*a3*b2**3*om*p3 + 3*a3**2*b2**3*om*p3 - 4*a4*b2**3*om*p3 - - 4*b2**2*b3*om*p3 + 14*a3*b2**2*b3*om*p3 + 15*b2*b3**2*om*p3 - - 28*b2**2*b4*om*p3 - 4*a3**3*b2**3*om**2*p3 + - 8*a3*a4*b2**3*om**2*p3 - 4*a5*b2**3*om**2*p3 - - a3**2*b2**2*b3*om**2*p3 - 2*a3*b2*b3**2*om**2*p3 - - b3**3*om**2*p3 + 4*a3*b2**2*b4*om**2*p3 + 4*b2*b3*b4*om**2*p3 - - 4*b2**2*b5*om**2*p3 - 80*b2**3*p4 + 4*b2**2*ll1*p4 - - 8*b2**3*om*p4 + 16*a3*b2**3*om*p4 - 16*b2**2*b3*om*p4 + - 3*a3**2*b2**3*om**2*p4 - 4*a4*b2**3*om**2*p4 + - 2*a3*b2**2*b3*om**2*p4 + 3*b2*b3**2*om**2*p4 - - 4*b2**2*b4*om**2*p4 + 40*b2**2*p0*xi + 8*b2**3*p0*xi - - 8*a3*b2**3*p0*xi - 8*a3**2*b2**3*p0*xi - 8*a3**3*b2**3*p0*xi + - 18*a3**4*b2**3*p0*xi + 16*a4*b2**3*p0*xi + 24*a3*a4*b2**3*p0*xi - - 76*a3**2*a4*b2**3*p0*xi + 40*a4**2*b2**3*p0*xi - - 24*a5*b2**3*p0*xi + 84*a3*a5*b2**3*p0*xi - 96*a6*b2**3*p0*xi - - 24*b2**2*b3*p0*xi + 8*a3*b2**2*b3*p0*xi + 8*a3**2*b2**2*b3*p0*xi - - 18*a3**3*b2**2*b3*p0*xi - 16*a4*b2**2*b3*p0*xi + - 58*a3*a4*b2**2*b3*p0*xi - 66*a5*b2**2*b3*p0*xi + - 40*b2**2*b4*p0*xi - 8*a3*b2**2*b4*p0*xi + - 18*a3**2*b2**2*b4*p0*xi - 40*a4*b2**2*b4*p0*xi - - 56*b2**2*b5*p0*xi - 18*a3*b2**2*b5*p0*xi - 24*b2**2*b6*p0*xi - - 32*b2**2*p1*xi - 16*b2**3*p1*xi + 8*a3*b2**3*p1*xi + - 8*a3**2*b2**3*p1*xi - 16*a3**3*b2**3*p1*xi - 16*a4*b2**3*p1*xi + - 52*a3*a4*b2**3*p1*xi - 60*a5*b2**3*p1*xi + 32*b2**2*b3*p1*xi - - 8*a3*b2**2*b3*p1*xi + 16*a3**2*b2**2*b3*p1*xi - - 36*a4*b2**2*b3*p1*xi - 48*b2**2*b4*p1*xi - 16*a3*b2**2*b4*p1*xi - - 20*b2**2*b5*p1*xi + 24*b2**2*p2*xi + 24*b2**3*p2*xi - - 8*a3*b2**3*p2*xi + 14*a3**2*b2**3*p2*xi - 32*a4*b2**3*p2*xi - - 40*b2**2*b3*p2*xi - 14*a3*b2**2*b3*p2*xi - 16*b2**2*b4*p2*xi - - 16*b2**2*p3*xi - 32*b2**3*p3*xi - 12*a3*b2**3*p3*xi - - 12*b2**2*b3*p3*xi + 8*b2**2*p4*xi - 8*b2**3*p4*xi)/(40.*b2**3*om) p = (s**alpha*(p0 + p1*s + p2*s**2 + p3*s**3 + p4*s**4 + p5*s**5))/ - dexp(om/s) p1p = (s**(-2.d0 + alpha)*(om*p0 + alpha*p0*s + om*p1*s + p1*s**2 + - alpha*p1*s**2 + om*p2*s**2 + 2*p2*s**3 + alpha*p2*s**3 + - om*p3*s**3 + 3*p3*s**4 + alpha*p3*s**4 + om*p4*s**4 + - 4*p4*s**5 + alpha*p4*s**5 + om*p5*s**5 + 5*p5*s**6 + - alpha*p5*s**6))/dexp(om/s) p2p = (s**(-4.d0 + alpha)*(om**2*p0 - 2*om*p0*s + 2*alpha*om*p0*s + - om**2*p1*s - alpha*p0*s**2 + alpha**2*p0*s**2 + - 2*alpha*om*p1*s**2 + om**2*p2*s**2 + alpha*p1*s**3 + - alpha**2*p1*s**3 + 2*om*p2*s**3 + 2*alpha*om*p2*s**3 + - om**2*p3*s**3 + 2*p2*s**4 + 3*alpha*p2*s**4 + alpha**2*p2*s**4 + - 4*om*p3*s**4 + 2*alpha*om*p3*s**4 + om**2*p4*s**4 + 6*p3*s**5 + - 5*alpha*p3*s**5 + alpha**2*p3*s**5 + 6*om*p4*s**5 + - 2*alpha*om*p4*s**5 + om**2*p5*s**5 + 12*p4*s**6 + - 7*alpha*p4*s**6 + alpha**2*p4*s**6 + 8*om*p5*s**6 + - 2*alpha*om*p5*s**6 + 20*p5*s**7 + 9*alpha*p5*s**7 + - alpha**2*p5*s**7))/dexp(om/s) c Print *,'alpha = ',alpha,' p0 = ',p0,' p1 = ',p1,' p2 = ',p2 c Print*,'p3 = ',p3,' p4 = ',p4,' p5 = ',p5 Return End