1      subroutine sp0001(gout)
2c $Id$
3c        *****  special fast routine for -p- loop for 0001 ****
4      implicit none
5      double precision gout(*)
6c
7      double precision a0, b0, c0, a1, b1, c1, a2
8      common/tabint/a0(333),b0(333),c0(334),a1(333),b1(333),c1(334),
9     +  a2(4000)
10c
11#include "auxvar.fh"
12#include "miscg.fh"
13#include "ginf.fh"
14#include "pgeom.fh"
15#include "shllfo.fh"
16#include "spgeom.fh"
17#include "qgeom.fh"
18#include "maxc.fh"
19c
20      double precision qq, theta, cq, dq
21      integer n
22      integer k, l, i
23      double precision sixty, tenm12, gout1, gout2, gout3, gout4
24      double precision q0, q1, ecd, qqq, qqtest
25      double precision dq00, dq01, dq10, dq11, v, aqx, aqz
26      double precision qperp, qperp2, cosp, sinp, p, q, pqab
27      double precision h0000, h0001, h0003, g0001, g0002, g0003
28      double precision theta2, theta3, theta4, u, g, t1, t2, t3
29c
30      data sixty,tenm12/60.0d0,1.0d-12/
31c
32      gout1 = 0.0d0
33      gout2 = 0.0d0
34      gout3 = 0.0d0
35      gout4 = 0.0d0
36      do 940 k = 1,ngc
37      gc = cg(k)
38      do 940 l = 1,ngd
39      gd = dg(l)
40      gcd = gc+gd
41      ecd = 1.0d0/gcd
42      cq = gd*ecd*rcd
43      dq = cq-rcd
44      qqq = cq*dq*gcd
45      if (qqq+sixty) 480,500,500
46  480 v = 0.0d0
47      go to 520
48  500 v =  dexp(qqq)*ecd
49  520 qqtest = cmaxc(k)*cmaxd(l)*v
50      if (qqtest-error1) 560,560,540
51  540 ismlq = 0
52      go to 600
53  560 if (qqtest-error2) 940,940,580
54  580 ismlq = 1
55  600 sc = csc(k)
56      sd = csd(l)
57      pc = cpc(k)
58      pd = cpd(l)
59      dq00 = sc*sd*v
60      dq01 = sc*pd*v
61      dq10 = pc*sd*v
62      dq11 = pc*pd*v
63      aqx = acx+sing*cq
64      aqz = acz+cosg*cq
65      qperp2 = aqx*aqx+acy2
66      qperp = dsqrt(qperp2)
67      if (qperp-tenm12) 640,640,620
68  620 cosp = -aqx/qperp
69      sinp = -acy/qperp
70      go to 660
71  640 cosp = 1.0d0
72      sinp = 0.0d0
73  660 h0000 = 0.d0
74      h0001 = 0.d0
75      h0003 = 0.d0
76      do 180 i = 1,ngangb
77      isml = ismlq+ismlp(i)
78      if (isml .ge. 2) go to 180
79      auxvar = var(isml+1)
80      pqab = aqz-app(i)
81      g = 1.d0/(ep(i)+ecd)
82      p = (pqab*pqab+qperp2)*g
83      if (p .le. auxvar) go to 140
84      q0 = dp00p(i)*dsqrt(0.7853981625d0/(p*(gp(i)+gcd)))
85      q1 = 0.5d0*q0/p
86      go to 160
87  140 q = dp00p(i)/dsqrt(gp(i)+gcd)
88      qq = p*12.5d0
89      n =  idint(qq)
90      theta = qq- dble(n)
91      theta2 = theta*(theta-1.d0)
92      theta3 = theta2*(theta-2.d0)
93      theta4 = theta2*(theta+1.d0)
94      q0 = (a0(n+1)+theta*b0(n+1)-theta3*c0(n+1)+theta4*c0(n+2))*q
95      q1 = (a1(n+1)+theta*b1(n+1)-theta3*c1(n+1)+theta4*c1(n+2))*q
96  160 u = g*q1
97      h0000 = h0000+q0
98      h0001 = h0001+u
99      h0003 = h0003-u*pqab
100  180 continue
101      h0001 = h0001*ecd*qperp
102      h0003 = h0003*ecd
103      p = dq*h0000
104      g0001 = h0001*cosp+p*sing
105      g0002 = h0001*sinp
106      g0003 = h0003+p*cosg
107      gout1 = gout1+dq00*h0000
108      gout2 = gout2+dq01*g0001
109      gout3 = gout3+dq01*g0002
110      gout4 = gout4+dq01*g0003
111 940  continue
112      gout(1) = gout1
113      t1 = gout2
114      t2 = gout3
115      t3 = gout4
116      gout(2) = p11*t1+p21*t2+p31*t3
117      gout(3) = p12*t1+p22*t2+p32*t3
118      gout(4) = p13*t1+p23*t2+p33*t3
119      return
120      end
121