ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/utils/interpolation.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/utils/interpolation.F90 (file contents):
Revision 2708 by gezelter, Fri Apr 14 19:57:04 2006 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 43 | Line 43
43   !!
44   !!  Created by Charles F. Vardeman II on 03 Apr 2006.
45   !!
46 < !!  PURPOSE: Generic Spline interplelation routines. These routines assume that we are on a uniform grid for
47 < !!           precomputation of spline parameters.
46 > !!  PURPOSE: Generic Spline interpolation routines.
47   !!
48   !! @author Charles F. Vardeman II
49 < !! @version $Id: interpolation.F90,v 1.1 2006-04-14 19:57:04 gezelter Exp $
49 > !! @version $Id: interpolation.F90,v 1.8 2006-05-17 15:37:15 gezelter Exp $
50  
51  
52 < module  INTERPOLATION
52 > module interpolation
53    use definitions
54    use status
55    implicit none
56    PRIVATE
57  
59 –  character(len = statusMsgSize) :: errMSG
60 –
58    type, public :: cubicSpline
59 <     private
60 <     integer :: np = 0
64 <     real(kind=dp) :: dx
59 >     logical :: isUniform = .false.
60 >     integer :: n = 0
61       real(kind=dp) :: dx_i
62       real (kind=dp), pointer,dimension(:)   :: x => null()
63 <     real (kind=dp), pointer,dimension(4,:) :: c => null()
63 >     real (kind=dp), pointer,dimension(:)   :: y => null()
64 >     real (kind=dp), pointer,dimension(:)   :: b => null()
65 >     real (kind=dp), pointer,dimension(:)   :: c => null()
66 >     real (kind=dp), pointer,dimension(:)   :: d => null()
67    end type cubicSpline
68  
69 <  interface splineLookup
71 <     module procedure multiSplint
72 <     module procedure splintd
73 <     module procedure splintd1
74 <     module procedure splintd2
75 <  end interface
76 <
77 <  interface newSpline
78 <     module procedure newSplineWithoutDerivs
79 <     module procedure newSplineWithDerivs
80 <  end interface
81 <
69 >  public :: newSpline
70    public :: deleteSpline
71 <
71 >  public :: lookupSpline
72 >  public :: lookupUniformSpline
73 >  public :: lookupNonuniformSpline
74 >  public :: lookupUniformSpline1d
75 >  
76   contains
77 +  
78  
79 <
80 <  subroutine newSplineWithoutDerivs(cs, x, y, yp1, ypn, boundary)
88 <
89 <    !************************************************************************
90 <    !
91 <    ! newSplineWithoutDerivs solves for slopes defining a cubic spline.
92 <    !
93 <    !  Discussion:
94 <    !
95 <    !    A tridiagonal linear system for the unknown slopes S(I) of
96 <    !    F at x(I), I=1,..., N, is generated and then solved by Gauss
97 <    !    elimination, with S(I) ending up in cs%C(2,I), for all I.
98 <    !
99 <    !  Reference:
100 <    !
101 <    !    Carl DeBoor,
102 <    !    A Practical Guide to Splines,
103 <    !    Springer Verlag.
104 <    !
105 <    !  Parameters:
106 <    !
107 <    !    Input, real x(N), the abscissas or X values of
108 <    !    the data points.  The entries of TAU are assumed to be
109 <    !    strictly increasing.
110 <    !
111 <    !    Input, real y(I), contains the function value at x(I) for
112 <    !      I = 1, N.
113 <    !
114 <    !    yp1 contains the slope at x(1) and ypn contains
115 <    !    the slope at x(N).
116 <    !
117 <    !    On output, the intermediate slopes at x(I) have been
118 <    !    stored in cs%C(2,I), for I = 2 to N-1.
119 <
79 >  subroutine newSpline(cs, x, y, isUniform)
80 >    
81      implicit none
82  
83      type (cubicSpline), intent(inout) :: cs
84      real( kind = DP ), intent(in) :: x(:), y(:)
85 <    real( kind = DP ), intent(in) :: yp1, ypn
86 <    character(len=*), intent(in) :: boundary
126 <    real( kind = DP ) :: g, divdif1, divdif3, dx
127 <    integer :: i, alloc_error, np
85 >    real( kind = DP ) :: fp1, fpn, p
86 >    REAL( KIND = DP), DIMENSION(size(x)-1) :: diff_y, H
87  
88 +    logical, intent(in) :: isUniform
89 +    integer :: i, alloc_error, n, k
90 +
91      alloc_error = 0
92  
93 <    if (cs%np .ne. 0) then
94 <       call handleWarning("interpolation::newSplineWithoutDerivs", &
95 <            "Type was already created")
93 >    if (cs%n .ne. 0) then
94 >       call handleWarning("interpolation::newSpline", &
95 >            "cubicSpline struct was already created")
96         call deleteSpline(cs)
97      end if
98  
99      ! make sure the sizes match
100  
101 <    if (size(x) .ne. size(y)) then
102 <       call handleError("interpolation::newSplineWithoutDerivs", &
101 >    n = size(x)
102 >    
103 >    if ( size(y) .ne. size(x) ) then
104 >       call handleError("interpolation::newSpline", &
105              "Array size mismatch")
106      end if
107 +    
108 +    cs%n = n
109 +    cs%isUniform = isUniform
110  
111 <    np = size(x)
145 <    cs%np = np
146 <
147 <    allocate(cs%x(np), stat=alloc_error)
111 >    allocate(cs%x(n), stat=alloc_error)
112      if(alloc_error .ne. 0) then
113 <       call handleError("interpolation::newSplineWithoutDerivs", &
113 >       call handleError("interpolation::newSpline", &
114              "Error in allocating storage for x")
115      endif
116  
117 <    allocate(cs%c(4,np), stat=alloc_error)
117 >    allocate(cs%y(n), stat=alloc_error)
118      if(alloc_error .ne. 0) then
119 <       call handleError("interpolation::newSplineWithoutDerivs", &
120 <            "Error in allocating storage for c")
119 >       call handleError("interpolation::newSpline", &
120 >            "Error in allocating storage for y")
121      endif
158 –      
159 –    do i = 1, np
160 –       cs%x(i) = x(i)
161 –       cs%c(1,i) = y(i)      
162 –    enddo
122  
123 <    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
124 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
125 <       cs%c(2,1) = yp1
126 <    else
168 <       cs%c(2,1) = 0.0_DP
123 >    allocate(cs%b(n), stat=alloc_error)
124 >    if(alloc_error .ne. 0) then
125 >       call handleError("interpolation::newSpline", &
126 >            "Error in allocating storage for b")
127      endif
128 <    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
129 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
130 <       cs%c(2,1) = ypn
131 <    else
132 <       cs%c(2,1) = 0.0_DP
128 >
129 >    allocate(cs%c(n), stat=alloc_error)
130 >    if(alloc_error .ne. 0) then
131 >       call handleError("interpolation::newSpline", &
132 >            "Error in allocating storage for c")
133      endif
134  
135 <    !
136 <    !  Set up the right hand side of the linear system.
137 <    !
138 <    do i = 2, cs%np - 1
139 <       cs%c(2,i) = 3.0_DP * ( &
182 <            (x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
183 <            (x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
184 <    end do
185 <    !
186 <    !  Set the diagonal coefficients.
187 <    !
188 <    cs%c(4,1) = 1.0_DP
189 <    do i = 2, cs%np - 1
190 <       cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
191 <    end do
192 <    cs%c(4,n) = 1.0_DP
193 <    !
194 <    !  Set the off-diagonal coefficients.
195 <    !
196 <    cs%c(3,1) = 0.0_DP
197 <    do i = 2, cs%np
198 <       cs%c(3,i) = x(i) - x(i-1)
199 <    end do
200 <    !
201 <    !  Forward elimination.
202 <    !
203 <    do i = 2, cs%np - 1
204 <       g = -cs%c(3,i+1) / cs%c(4,i-1)
205 <       cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
206 <       cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
207 <    end do
208 <    !
209 <    !  Back substitution for the interior slopes.
210 <    !
211 <    do i = cs%np - 1, 2, -1
212 <       cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
213 <    end do
214 <    !
215 <    !  Now compute the quadratic and cubic coefficients used in the
216 <    !  piecewise polynomial representation.
217 <    !
218 <    do i = 1, cs%np - 1
219 <       dx = x(i+1) - x(i)
220 <       divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
221 <       divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
222 <       cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
223 <       cs%c(4,i) = divdif3 / ( dx * dx )
224 <    end do
135 >    allocate(cs%d(n), stat=alloc_error)
136 >    if(alloc_error .ne. 0) then
137 >       call handleError("interpolation::newSpline", &
138 >            "Error in allocating storage for d")
139 >    endif
140  
141 <    cs%c(3,np) = 0.0_DP
227 <    cs%c(4,np) = 0.0_DP
141 >    ! make sure we are monotinically increasing in x:
142  
143 <    cs%dx = dx
144 <    cs%dxi = 1.0_DP / dx
145 <    return
146 <  end subroutine newSplineWithoutDerivs
143 >    h = diff(x)
144 >    if (any(h <= 0)) then
145 >       call handleError("interpolation::newSpline", &
146 >            "Negative dx interval found")
147 >    end if
148  
149 <  subroutine newSplineWithDerivs(cs, x, y, yp)
149 >    ! load x and y values into the cubicSpline structure:
150  
151 <    !************************************************************************
152 <    !
153 <    ! newSplineWithDerivs
151 >    do i = 1, n
152 >       cs%x(i) = x(i)
153 >       cs%y(i) = y(i)
154 >    end do
155  
156 <    implicit none
156 >    ! Calculate coefficients for the tridiagonal system: store
157 >    ! sub-diagonal in B, diagonal in D, difference quotient in C.
158  
159 <    type (cubicSpline), intent(inout) :: cs
160 <    real( kind = DP ), intent(in) :: x(:), y(:), yp(:)
161 <    real( kind = DP ) :: g, divdif1, divdif3, dx
245 <    integer :: i, alloc_error, np
159 >    cs%b(1:n-1) = h
160 >    diff_y = diff(y)
161 >    cs%c(1:n-1) = diff_y / h
162  
163 <    alloc_error = 0
164 <
165 <    if (cs%np .ne. 0) then
166 <       call handleWarning("interpolation::newSplineWithDerivs", &
167 <            "Type was already created")
168 <       call deleteSpline(cs)
163 >    if (n == 2) then
164 >       ! Assume the derivatives at both endpoints are zero
165 >       ! another assumption could be made to have a linear interpolant
166 >       ! between the two points.  In that case, the b coefficients
167 >       ! below would be diff_y(1)/h(1) and the c and d coefficients would
168 >       ! both be zero.
169 >       cs%b(1) = 0.0_dp
170 >       cs%c(1) = -3.0_dp * (diff_y(1)/h(1))**2
171 >       cs%d(1) = -2.0_dp * (diff_y(1)/h(1))**3
172 >       cs%b(2) = cs%b(1)
173 >       cs%c(2) = 0.0_dp
174 >       cs%d(2) = 0.0_dp
175 >       cs%dx_i = 1.0_dp / h(1)
176 >      return
177      end if
178  
179 <    ! make sure the sizes match
180 <
181 <    if ((size(x) .ne. size(y)).or.(size(x) .ne. size(yp))) then
182 <       call handleError("interpolation::newSplineWithDerivs", &
183 <            "Array size mismatch")
260 <    end if
261 <    
262 <    np = size(x)
263 <    cs%np = np
179 >    cs%d(1) = 2.0_dp * cs%b(1)
180 >    do i = 2, n-1
181 >      cs%d(i) = 2.0_dp * (cs%b(i) + cs%b(i-1))
182 >    end do
183 >    cs%d(n) = 2.0_dp * cs%b(n-1)
184  
185 <    allocate(cs%x(np), stat=alloc_error)
186 <    if(alloc_error .ne. 0) then
267 <       call handleError("interpolation::newSplineWithDerivs", &
268 <            "Error in allocating storage for x")
269 <    endif
185 >    ! Calculate estimates for the end slopes using polynomials
186 >    ! that interpolate the data nearest the end.
187      
188 <    allocate(cs%c(4,np), stat=alloc_error)
189 <    if(alloc_error .ne. 0) then
190 <       call handleError("interpolation::newSplineWithDerivs", &
191 <            "Error in allocating storage for c")
192 <    endif
193 <    
194 <    do i = 1, np
195 <       cs%x(i) = x(i)
196 <       cs%c(1,i) = y(i)      
197 <       cs%c(2,i) = yp(i)
198 <    enddo
199 <    !
200 <    !  Set the diagonal coefficients.
201 <    !
202 <    cs%c(4,1) = 1.0_DP
203 <    do i = 2, cs%np - 1
204 <       cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
188 >    fp1 = cs%c(1) - cs%b(1)*(cs%c(2) - cs%c(1))/(cs%b(1) + cs%b(2))
189 >    if (n > 3) then
190 >      fp1 = fp1 + cs%b(1)*((cs%b(1) + cs%b(2))*(cs%c(3) - cs%c(2))/ &
191 >           (cs%b(2) + cs%b(3)) - cs%c(2) + cs%c(1))/(x(4) - x(1))
192 >    end if
193 >          
194 >    fpn = cs%c(n-1) + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2))/(cs%b(n-2) + cs%b(n-1))
195 >    if (n > 3) then
196 >      fpn = fpn + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2) - (cs%b(n-2) + cs%b(n-1))* &
197 >           (cs%c(n-2) - cs%c(n-3))/(cs%b(n-2) + cs%b(n-3)))/(x(n) - x(n-3))
198 >    end if
199 >
200 >    ! Calculate the right hand side and store it in C.
201 >
202 >    cs%c(n) = 3.0_dp * (fpn - cs%c(n-1))
203 >    do i = n-1,2,-1
204 >      cs%c(i) = 3.0_dp * (cs%c(i) - cs%c(i-1))
205      end do
206 <    cs%c(4,n) = 1.0_DP
207 <    !
208 <    !  Set the off-diagonal coefficients.
209 <    !
210 <    cs%c(3,1) = 0.0_DP
211 <    do i = 2, cs%np
212 <       cs%c(3,i) = x(i) - x(i-1)
206 >    cs%c(1) = 3.0_dp * (cs%c(1) - fp1)
207 >
208 >    ! Solve the tridiagonal system.
209 >
210 >    do k = 2, n
211 >      p = cs%b(k-1) / cs%d(k-1)
212 >      cs%d(k) = cs%d(k) - p*cs%b(k-1)
213 >      cs%c(k) = cs%c(k) - p*cs%c(k-1)
214      end do
215 <    !
216 <    !  Forward elimination.
217 <    !
300 <    do i = 2, cs%np - 1
301 <       g = -cs%c(3,i+1) / cs%c(4,i-1)
302 <       cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
303 <       cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
215 >    cs%c(n) = cs%c(n) / cs%d(n)
216 >    do k = n-1, 1, -1
217 >      cs%c(k) = (cs%c(k) - cs%b(k) * cs%c(k+1)) / cs%d(k)
218      end do
305 –    !
306 –    !  Back substitution for the interior slopes.
307 –    !
308 –    do i = cs%np - 1, 2, -1
309 –       cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
310 –    end do
311 –    !
312 –    !  Now compute the quadratic and cubic coefficients used in the
313 –    !  piecewise polynomial representation.
314 –    !
315 –    do i = 1, cs%np - 1
316 –       dx = x(i+1) - x(i)
317 –       divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
318 –       divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
319 –       cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
320 –       cs%c(4,i) = divdif3 / ( dx * dx )
321 –    end do
219  
220 <    cs%c(3,np) = 0.0_DP
324 <    cs%c(4,np) = 0.0_DP
220 >    ! Calculate the coefficients defining the spline.
221  
222 <    cs%dx = dx
223 <    cs%dxi = 1.0_DP / dx
222 >    cs%d(1:n-1) = diff(cs%c) / (3.0_dp * h)
223 >    cs%b(1:n-1) = diff_y / h - h * (cs%c(1:n-1) + h * cs%d(1:n-1))
224 >    cs%b(n) = cs%b(n-1) + h(n-1) * (2.0_dp*cs%c(n-1) + h(n-1)*3.0_dp*cs%d(n-1))
225  
226 <    return
227 <  end subroutine newSplineWithoutDerivs
226 >    if (isUniform) then
227 >       cs%dx_i = 1.0_dp / (x(2) - x(1))
228 >    endif
229  
230 +    return
231 +    
232 +  contains
233 +    
234 +    function diff(v)
235 +      ! Auxiliary function to compute the forward difference
236 +      ! of data stored in a vector v.
237 +      
238 +      implicit none
239 +      real (kind = dp), dimension(:), intent(in) :: v
240 +      real (kind = dp), dimension(size(v)-1) :: diff
241 +      
242 +      integer :: n
243 +      
244 +      n = size(v)
245 +      diff = v(2:n) - v(1:n-1)
246 +      return
247 +    end function diff
248 +    
249 +  end subroutine newSpline
250 +      
251    subroutine deleteSpline(this)
252  
253      type(cubicSpline) :: this
# Line 342 | Line 261 | contains
261         this%c => null()
262      end if
263      
264 <    this%np = 0
264 >    this%n = 0
265      
266    end subroutine deleteSpline
267  
268 <  subroutine lookup_nonuniform_spline(cs, xval, yval)
268 >  subroutine lookupNonuniformSpline(cs, xval, yval)
269      
351 –    !*************************************************************************
352 –    !
353 –    ! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant.
354 –    !
355 –    !  Discussion:
356 –    !
357 –    !    newSpline must be called first, to set up the
358 –    !    spline data from the raw function and derivative data.
359 –    !
360 –    !  Modified:
361 –    !
362 –    !    06 April 1999
363 –    !
364 –    !  Reference:
365 –    !
366 –    !    Conte and de Boor,
367 –    !    Algorithm PCUBIC,
368 –    !    Elementary Numerical Analysis,
369 –    !    1973, page 234.
370 –    !
371 –    !  Parameters:
372 –    !
270      implicit none
271  
272      type (cubicSpline), intent(in) :: cs
273      real( kind = DP ), intent(in)  :: xval
274      real( kind = DP ), intent(out) :: yval
275 +    real( kind = DP ) ::  dx
276      integer :: i, j
277      !
278      !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
279      !  or is nearest to xval.
280      !
281 <    j = cs%np - 1
281 >    j = cs%n - 1
282  
283 <    do i = 1, cs%np - 2
283 >    do i = 0, cs%n - 2
284  
285         if ( xval < cs%x(i+1) ) then
286            j = i
# Line 394 | Line 292 | contains
292      !  Evaluate the cubic polynomial.
293      !
294      dx = xval - cs%x(j)
295 <
398 <    yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
295 >    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
296      
297      return
298 <  end subroutine lookup_nonuniform_spline
298 >  end subroutine lookupNonuniformSpline
299  
300 <  subroutine lookup_uniform_spline(cs, xval, yval)
300 >  subroutine lookupUniformSpline(cs, xval, yval)
301      
405 –    !*************************************************************************
406 –    !
407 –    ! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant.
408 –    !
409 –    !  Discussion:
410 –    !
411 –    !    newSpline must be called first, to set up the
412 –    !    spline data from the raw function and derivative data.
413 –    !
414 –    !  Modified:
415 –    !
416 –    !    06 April 1999
417 –    !
418 –    !  Reference:
419 –    !
420 –    !    Conte and de Boor,
421 –    !    Algorithm PCUBIC,
422 –    !    Elementary Numerical Analysis,
423 –    !    1973, page 234.
424 –    !
425 –    !  Parameters:
426 –    !
302      implicit none
303  
304      type (cubicSpline), intent(in) :: cs
305      real( kind = DP ), intent(in)  :: xval
306      real( kind = DP ), intent(out) :: yval
307 +    real( kind = DP ) ::  dx
308      integer :: i, j
309      !
310      !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
311      !  or is nearest to xval.
312 +    
313 +    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
314 +    
315 +    dx = xval - cs%x(j)
316 +    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
317 +    
318 +    return
319 +  end subroutine lookupUniformSpline
320  
321 <    j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dxi) + 1))
321 >  subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
322 >    
323 >    implicit none
324  
325 +    type (cubicSpline), intent(in) :: cs
326 +    real( kind = DP ), intent(in)  :: xval
327 +    real( kind = DP ), intent(out) :: yval, dydx
328 +    real( kind = DP ) :: dx
329 +    integer :: i, j
330 +    
331 +    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
332 +    !  or is nearest to xval.
333 +
334 +
335 +    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
336 +    
337      dx = xval - cs%x(j)
338 +    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
339  
340 <    yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
340 >    dydx = cs%b(j) + dx*(2.0_dp * cs%c(j) + 3.0_dp * dx * cs%d(j))
341 >      
342 >    return
343 >  end subroutine lookupUniformSpline1d
344 >
345 >  subroutine lookupSpline(cs, xval, yval)
346 >
347 >    type (cubicSpline), intent(in) :: cs
348 >    real( kind = DP ), intent(inout) :: xval
349 >    real( kind = DP ), intent(inout) :: yval
350      
351 +    if (cs%isUniform) then
352 +       call lookupUniformSpline(cs, xval, yval)
353 +    else
354 +       call lookupNonuniformSpline(cs, xval, yval)
355 +    endif
356 +
357      return
358 <  end subroutine lookup_uniform_spline
358 >  end subroutine lookupSpline
359    
360 < end module INTERPOLATION
360 > end module interpolation

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines