← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_242195/PROC_242195/Source/PDF/jeppe02.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_242195/PROC_242195/Source/PDF/jeppe02.f
 
1
      subroutine jeppe1(nx,my,xx,yy,ff,cc)
 
2
      implicit real*8(a-h,o-z)
 
3
      parameter(nnx=49,mmy=37)
 
4
      dimension xx(nx),yy(my),ff(nnx,mmy),ff1(nnx,mmy),ff2(nnx,mmy),
 
5
     xff12(nnx,mmy),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
 
6
     xcl(16),cc(nx,my,4,4),iwt(16,16)
 
7
 
 
8
      data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 
9
     x            0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
 
10
     x            -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
 
11
     x            2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
 
12
     x            0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
 
13
     x            0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
 
14
     x            0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
 
15
     x            0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
 
16
     x            -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
 
17
     x            0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
 
18
     x            9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
 
19
     x            -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
 
20
     x            2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
 
21
     x            0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
 
22
     x            -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
 
23
     x            4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
 
24
 
 
25
 
 
26
      do 42 m=1,my
 
27
      dx=xx(2)-xx(1)
 
28
      ff1(1,m)=(ff(2,m)-ff(1,m))/dx
 
29
      dx=xx(nx)-xx(nx-1)
 
30
      ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
 
31
      do 41 n=2,nx-1
 
32
      ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
 
33
     xff(n+1,m))
 
34
   41 continue
 
35
   42 continue
 
36
 
 
37
      do 44 n=1,nx
 
38
      dy=yy(2)-yy(1)
 
39
      ff2(n,1)=(ff(n,2)-ff(n,1))/dy
 
40
      dy=yy(my)-yy(my-1)
 
41
      ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
 
42
      do 43 m=2,my-1
 
43
      ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
 
44
     xff(n,m+1))
 
45
   43 continue
 
46
   44 continue
 
47
 
 
48
      do 46 m=1,my
 
49
      dx=xx(2)-xx(1)
 
50
      ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
 
51
      dx=xx(nx)-xx(nx-1)
 
52
      ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
 
53
      do 45 n=2,nx-1
 
54
      ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
 
55
     xff2(n+1,m))
 
56
   45 continue
 
57
   46 continue
 
58
 
 
59
      do 53 n=1,nx-1
 
60
      do 52 m=1,my-1
 
61
      d1=xx(n+1)-xx(n)
 
62
      d2=yy(m+1)-yy(m)
 
63
      d1d2=d1*d2
 
64
 
 
65
      yy0(1)=ff(n,m)
 
66
      yy0(2)=ff(n+1,m)
 
67
      yy0(3)=ff(n+1,m+1)
 
68
      yy0(4)=ff(n,m+1)
 
69
 
 
70
      yy1(1)=ff1(n,m)
 
71
      yy1(2)=ff1(n+1,m)
 
72
      yy1(3)=ff1(n+1,m+1)
 
73
      yy1(4)=ff1(n,m+1)
 
74
 
 
75
      yy2(1)=ff2(n,m)
 
76
      yy2(2)=ff2(n+1,m)
 
77
      yy2(3)=ff2(n+1,m+1)
 
78
      yy2(4)=ff2(n,m+1)
 
79
 
 
80
      yy12(1)=ff12(n,m)
 
81
      yy12(2)=ff12(n+1,m)
 
82
      yy12(3)=ff12(n+1,m+1)
 
83
      yy12(4)=ff12(n,m+1)
 
84
 
 
85
      do 47 k=1,4
 
86
      z(k)=yy0(k)
 
87
      z(k+4)=yy1(k)*d1
 
88
      z(k+8)=yy2(k)*d2
 
89
      z(k+12)=yy12(k)*d1d2
 
90
   47 continue
 
91
 
 
92
      do 49 l=1,16
 
93
      xxd=0.
 
94
      do 48 k=1,16
 
95
      xxd=xxd+iwt(k,l)*z(k)
 
96
   48 continue
 
97
      cl(l)=xxd
 
98
   49 continue
 
99
      l=0
 
100
      do 51 k=1,4
 
101
      do 50 j=1,4
 
102
      l=l+1
 
103
      cc(n,m,k,j)=cl(l)
 
104
   50 continue
 
105
   51 continue
 
106
   52 continue
 
107
   53 continue
 
108
      return
 
109
      end
 
110
 
 
111
      subroutine jeppe2(x,y,nx,my,xx,yy,cc,z)
 
112
      implicit real*8(a-h,o-z)
 
113
      dimension xx(nx),yy(my),cc(nx,my,4,4)      
 
114
 
 
115
      n=locx(xx,nx,x)
 
116
      m=locx(yy,my,y)
 
117
 
 
118
      t=(x-xx(n))/(xx(n+1)-xx(n))
 
119
      u=(y-yy(m))/(yy(m+1)-yy(m))
 
120
 
 
121
      z=0.
 
122
      do 1 l=4,1,-1
 
123
      z=t*z+((cc(n,m,l,4)*u+cc(n,m,l,3))*u
 
124
     .       +cc(n,m,l,2))*u+cc(n,m,l,1)
 
125
    1 continue
 
126
      return
 
127
      end
 
128
 
 
129
      integer function locx(xx,nx,x)
 
130
      implicit real*8(a-h,o-z)
 
131
      dimension xx(nx)
 
132
      if(x.le.xx(1)) then
 
133
      locx=1
 
134
      return
 
135
      endif
 
136
      if(x.ge.xx(nx)) then 
 
137
      locx=nx-1  
 
138
      return
 
139
      endif
 
140
      ju=nx+1
 
141
      jl=0
 
142
    1 if((ju-jl).le.1) go to 2
 
143
      jm=(ju+jl)/2
 
144
      if(x.ge.xx(jm)) then
 
145
      jl=jm
 
146
      else
 
147
      ju=jm
 
148
      endif
 
149
      go to 1
 
150
    2 locx=jl
 
151
      return
 
152
      end
 
153
 
 
154
 
 
155
      real*8 function  polderiv(x1,x2,x3,y1,y2,y3)
 
156
      implicit real*8(a-h,o-z)
 
157
      polderiv=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
 
158
     .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
 
159
      return
 
160
      end