1C Output from Public domain Ratfor, version 1.05
2      subroutine sakj(x,z,p,iker,dens,psi,score,nx,nz,h,alpha,kappa,xlam
3     *)
4      double precision dens(nz),score(nz),psi(nz),h,kappa
5      double precision z(nz),x(nx),xlam(nx),p(nx),qrange,pi
6      double precision con1,sum,sqsum,xsd,a,fifth,hinv,half
7      double precision xn,xker,dxker,ddxker,fact,xponen,alpha,glog,zero,
8     *one,two
9      parameter( zero = 0.d0)
10      parameter( one = 1.d0)
11      parameter( two = 2.d0)
12      parameter( four = 4.d0)
13      parameter( half = 0.5d0)
14      parameter( fifth = 0.2d0)
15      parameter( pi = 3.141593d0)
16      xn=nx
17      if(iker.eq.0)then
18      con1=one/sqrt(2.0*pi)
19      else
20      if(iker.eq.1)then
21      con1=one/pi
22      endif
23      endif
24      if(h.le.0.)then
25      sum=0.
26      sqsum=0.
27      do23006 i=1,nx
28      sqsum=sqsum+x(i)*x(i)*p(i)
29      sum=sum+x(i)*p(i)
3023006 continue
3123007 continue
32      xsd=dsqrt(sqsum-sum*sum)
33      sum=zero
34      i=1
3523008 if(.not.(i.lt.nx))goto 23010
36      sum=sum+p(i)
37      if(sum.lt..25)then
38      goto 23009
39      else
40      qrange=x(i)
41      goto 23010
42      endif
4323009 i=i+1
44      goto 23008
4523010 continue
46      sum=one
47      i=nx
4823013 if(.not.(i.gt.0))goto 23015
49      sum=sum-p(i)
50      if(sum.gt..75)then
51      goto 23014
52      else
53      qrange=x(i)-qrange
54      goto 23015
55      endif
5623014 i=i-1
57      goto 23013
5823015 continue
59      a=min(xsd,qrange/1.34)
60      h=kappa*a/(xn**fifth)
61      endif
62      hinv=one/h
63      do23018 j=1,nx
64      xker=0.
65      if(iker.eq.0)then
66      do23022 i=1,nx
67      xponen=(x(j)-x(i))*hinv
68      xponen=half*xponen**2
69      xker=xker+p(i)*exp(-xponen)*hinv
7023022 continue
7123023 continue
72      else
73      if(iker.eq.1)then
74      do23026 i=1,nx
75      xponen=(x(j)-x(i))*hinv
76      xker=xker+p(i)*hinv/(1+xponen**2)
7723026 continue
7823027 continue
79      endif
80      endif
81      xlam(j)=con1*xker
8223018 continue
8323019 continue
84      glog=zero
85      do23028 i=1,nx
86      glog=glog+p(i)*log(xlam(i))
8723028 continue
8823029 continue
89      g=exp(glog)
90      ginv=one/g
91      do23030 i=1,nx
92      xlam(i)=hinv/((xlam(i)*ginv)**(-alpha))
9323030 continue
9423031 continue
95      do23032 j=1,nz
96      xker=zero
97      dxker=zero
98      ddxker=zero
99      if(iker.eq.0)then
100      do23036 i=1,nx
101      xponen=(z(j)-x(i))*xlam(i)
102      fact=exp(-half*xponen*xponen)*xlam(i)
103      xker=xker+p(i)*fact
104      dxker=dxker-p(i)*fact*xponen*xlam(i)
105      ddxker=ddxker- p(i)*fact*(one - xponen**2)*xlam(i)**2
10623036 continue
10723037 continue
108      else
109      if(iker.eq.1)then
110      do23040 i=1,nx
111      xponen=(z(j)-x(i))*xlam(i)
112      fact=xlam(i)/(one+xponen**2)
113      xker=xker+p(i)*fact
114      dxker=dxker-p(i)*two*xponen*fact**2
115      ddxker=ddxker- p(i)*two*(fact**2)*(xlam(i)- four*(xponen**2)*fact)
11623040 continue
11723041 continue
118      endif
119      endif
120      dens(j)=con1*xker
121      psi(j)=-(dxker/xker)
122      score(j)=(dxker/xker)**2-ddxker/xker
12323032 continue
12423033 continue
125      return
126      end
127