← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_141512/PROC_141512/Source/getissud.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_141512/PROC_141512/Source/getissud.f
 
1
C...GETISSUD performs an interpolation/extrapolation in 3 dimensions by
 
2
C...fitting quadratic splines using 4 points in each dimension
 
3
      double precision function getissud(ibeam,kfl,x1,x2,pt2)
 
4
      implicit none
 
5
 
 
6
      include 'sudgrid.inc'
 
7
 
 
8
c   Arguments
 
9
      integer ibeam,kfl
 
10
      double precision x1,x2,pt2
 
11
c   Storing values for the interpolation
 
12
      double precision smallgrid(4,4),minigrid(4) ! pt2,x1
 
13
c   Local variables
 
14
      integer ipt2,ix1,ix2,ilo,ihi,i,j,k,kkfl,ipoints
 
15
      double precision pt2i,x2i,x1i,minpoint,maxpoint,x(3)
 
16
      integer nerr
 
17
      data nerr/0/
 
18
 
 
19
      getissud=0
 
20
 
 
21
      x(1)=log(x2)
 
22
      x(2)=x1
 
23
      x(3)=log(pt2)
 
24
 
 
25
      kkfl=kfl
 
26
      if(ibeam.lt.0) kkfl=-kkfl
 
27
      if(kkfl.lt.-2) kkfl=iabs(kfl)
 
28
      if(iabs(kkfl).eq.21) kkfl=0
 
29
      if(kkfl.eq.5) then
 
30
        ipoints=2
 
31
      else
 
32
        ipoints=1
 
33
      endif        
 
34
      if(kkfl.gt.5) then
 
35
        if(nerr.lt.10)
 
36
     $     write(*,*)'GETISSUD Warning: flavor ',kfl,' not supported'
 
37
        nerr=nerr+1
 
38
        getissud=1
 
39
        return
 
40
      endif
 
41
 
 
42
      if(x(1).lt.points(1,ipoints).or.
 
43
     $   x(1).gt.points(nx2,ipoints).and.nerr.lt.10)
 
44
     $   then
 
45
        write(*,*) 'GETISSUD Warning: extrapolation in x2: ',x2
 
46
        nerr=nerr+1
 
47
      endif
 
48
 
 
49
      if(x(2).lt.points(nx2+1,ipoints).or.
 
50
     $   x(2).gt.points(nx2+nx1,ipoints)
 
51
     $   .and.nerr.lt.10) then
 
52
        write(*,*) 'GETISSUD Warning: extrapolation in x1: ',x1
 
53
        nerr=nerr+1
 
54
      endif
 
55
 
 
56
      if(kkfl.eq.5.and.pt2.lt.22.3109)then
 
57
        getissud=1d0
 
58
        return
 
59
      endif
 
60
 
 
61
      if(kkfl.eq.5.and.x1.gt.0.6)then
 
62
        getissud=0d0
 
63
        return
 
64
      endif
 
65
 
 
66
      if(x(3).lt.points(nx2+nx1+1,ipoints)) then
 
67
        write(*,*) 'GETISSUD Error! pt2 = ',exp(x(3)),' < ',
 
68
     $     exp(points(nx2+nx1+1,ipoints)),' = min(pt2). Not allowed!'
 
69
        write(*,*) 'You need to regenerate grid with new pt2min.'
 
70
        stop
 
71
      endif
 
72
 
 
73
      if(x(3).lt.points(nx2+nx1+1,ipoints).or.
 
74
     $   x(3).gt.points(nx2+nx1+npt2,ipoints)
 
75
     $   .and.nerr.lt.10) then
 
76
        write(*,*) 'GETISSUD Warning: extrapolation in pt2: ',pt2
 
77
        nerr=nerr+1
 
78
      endif
 
79
 
 
80
 
 
81
c   Find nearest points by binary method
 
82
c   x2
 
83
      ilo=1
 
84
      ihi=nx2
 
85
      do while(ihi.gt.ilo+1)
 
86
        ix2=ilo+(ihi-ilo)/2
 
87
        if(x(1).gt.points(ix2,ipoints))then
 
88
          ilo=ix2
 
89
        else
 
90
          ihi=ix2
 
91
        endif
 
92
      enddo
 
93
      if(x(1).lt.points(ix2,ipoints))
 
94
     $   ix2=ix2-1
 
95
      ix2=max(2,min(ix2,nx2-2))
 
96
 
 
97
c      print *,'x2: ',ix2,x(1),(points(i,ipoints),i=ix2-1,ix2+2)
 
98
 
 
99
c   x1
 
100
      ilo=1
 
101
      ihi=nx1
 
102
      do while(ihi.gt.ilo+1)
 
103
        ix1=ilo+(ihi-ilo)/2
 
104
        if(x(2).gt.points(nx2+ix1,ipoints))then
 
105
          ilo=ix1
 
106
        else
 
107
          ihi=ix1
 
108
        endif
 
109
      enddo
 
110
      if(x(2).lt.points(nx2+ix1,ipoints))
 
111
     $   ix1=ix1-1
 
112
      ix1=max(2,min(ix1,nx1-2))
 
113
 
 
114
      do while(kkfl.eq.5.and.
 
115
     $   points(nx2+ix1+2,ipoints).gt.0.6)
 
116
        ix1=ix1-1
 
117
      enddo
 
118
 
 
119
c      print *,'x1: ',ix1,x(2),(points(nx2+i,ipoints),i=ix1-1,ix1+2)
 
120
 
 
121
c   pt2
 
122
      ilo=1
 
123
      ihi=npt2
 
124
      do while(ihi.gt.ilo+1)
 
125
        ipt2=ilo+(ihi-ilo)/2
 
126
        if(x(3).gt.points(nx2+nx1+ipt2,ipoints))then
 
127
          ilo=ipt2
 
128
        else
 
129
          ihi=ipt2
 
130
        endif
 
131
      enddo
 
132
      if(x(3).lt.points(nx2+nx1+ipt2,ipoints))
 
133
     $   ipt2=ipt2-1
 
134
      ipt2=max(2,min(ipt2,npt2-2))
 
135
 
 
136
      do while(kkfl.eq.5.and.
 
137
     $   exp(points(nx2+nx1+ipt2-1,ipoints)).lt.22.3109)
 
138
        ipt2=ipt2+1
 
139
      enddo
 
140
c      print *,'pt2: ',ipt2,x(3),(points(nx2+nx1+i,ipoints),i=ipt2-1,ipt2+2)
 
141
c      print *,'pt: ',ipt2,exp(x(3)/2),
 
142
c     $   (exp(points(nx2+nx1+i,ipoints)/2),i=ipt2-1,ipt2+2)
 
143
 
 
144
C   Now perform inter-/extra-polation
 
145
 
 
146
C   Start with x2, which should have the flattest distribution
 
147
C   Calculate sud(x2,ax1,apt2) for the 4x4 apt2 and ax1
 
148
C   Then continue with pt2 and calculate sud(x2,ax1,pt2)
 
149
C   for the 4 ax1
 
150
C   Finally calculate sud(x2,x1,pt2)
 
151
 
 
152
      do i=1,4
 
153
        do j=1,4
 
154
c          print *,'x1,pt:',points(nx2+ix1-2+i,ipoints),
 
155
c     $       exp(points(nx2+nx1+ipt2-2+j,ipoints)/2)
 
156
          call splint2(sudgrid(ix2-1,ix1-2+i,ipt2-2+j,kkfl),
 
157
     $       points(ix2-1,ipoints),4,x(1),smallgrid(j,i))
 
158
          smallgrid(j,i)=max(0d0,min(1d0,smallgrid(j,i)))
 
159
        enddo
 
160
      enddo
 
161
 
 
162
      do i=1,4
 
163
        call splint2(smallgrid(1,i),
 
164
     $     points(nx2+nx1+ipt2-1,ipoints),4,x(3),minigrid(i))
 
165
          minigrid(i)=max(0d0,min(1d0,minigrid(i)))
 
166
      enddo
 
167
 
 
168
      call splint2(minigrid,
 
169
     $   points(nx2+ix1-1,ipoints),4,x(2),getissud)
 
170
      getissud=max(0d0,min(1d0,getissud))
 
171
      
 
172
c      print *,'Result: ',getissud
 
173
 
 
174
      return
 
175
      end
 
176
 
 
177
 
 
178
      subroutine splint2(ypoints,xpoints,npoints,x,ans)
 
179
      implicit none
 
180
 
 
181
C   arguments
 
182
      integer npoints
 
183
      double precision ypoints(npoints),xpoints(npoints)
 
184
      double precision x,ans
 
185
C   local variables
 
186
      double precision a0,a1,a2,sd
 
187
      integer ifail,i
 
188
 
 
189
      CALL DLSQP2(npoints,xpoints,ypoints,a0,a1,a2,sd,ifail)
 
190
 
 
191
c      print *,'Point, interpolation:'
 
192
c      do i=1,npoints
 
193
c        print *,exp(xpoints(i)),ypoints(i),
 
194
c     $     a0+a1*xpoints(i)+a2*xpoints(i)**2
 
195
c      enddo
 
196
 
 
197
      ans=a0+a1*x+a2*x**2
 
198
c      print *,x,ans
 
199
 
 
200
      return
 
201
      end