@@ -78,7 +78,18 @@ module fitpack_surfaces
7878 procedure :: new_fit = > surf_new_fit
7979
8080 ! > Generate/update fitting curve, with optional smoothing
81- procedure :: fit = > surface_fit_automatic_knots
81+ procedure :: fit = > surface_fit_automatic_knots
82+ procedure :: interpolate = > surface_fit_interpolating
83+ procedure :: least_squares = > surface_fit_least_squares
84+
85+ ! > Evaluate surface at given coordinates
86+ procedure , private :: surface_eval_one
87+ procedure , private :: surface_eval_many
88+ generic :: eval = > surface_eval_one,surface_eval_many
89+
90+ ! > Evaluate surface on a grid
91+ procedure , private :: surface_eval_gridded
92+ generic :: eval_ongrid = > surface_eval_gridded
8293
8394 ! > Evaluate derivatives at given coordinates
8495 procedure , private :: surface_derivatives_gridded
@@ -87,6 +98,9 @@ module fitpack_surfaces
8798 generic :: dfdx = > surface_derivatives_one,surface_derivatives_many
8899 generic :: dfdx_ongrid = > surface_derivatives_gridded
89100
101+ ! > Properties: MSE
102+ procedure , non_overridable :: mse = > surface_error
103+
90104 end type fitpack_surface
91105
92106 interface fitpack_surface
@@ -144,6 +158,104 @@ integer(FP_FLAG) function surface_fit_automatic_knots(this,smoothing,order) resu
144158 end function surface_fit_automatic_knots
145159
146160
161+ ! Find interpolating surface
162+ integer (FP_FLAG) function surface_fit_interpolating(this) result(ierr)
163+ class(fitpack_surface), intent (inout ) :: this
164+ ierr = surface_fit_automatic_knots(this,smoothing= zero)
165+ end function surface_fit_interpolating
166+
167+ ! Fit a surface to least squares of the current knots
168+ integer (FP_FLAG) function surface_fit_least_squares(this) result(ierr)
169+ class(fitpack_surface), intent (inout ) :: this
170+ this% iopt = IOPT_NEW_LEASTSQUARES
171+ ierr = this% fit()
172+ end function surface_fit_least_squares
173+
174+ ! > Evaluate surface on a list of (x(i),y(i)) scattered points using bispeu
175+ function surface_eval_many (this ,x ,y ,ierr ) result(f)
176+ class(fitpack_surface), intent (inout ) :: this
177+ real (FP_REAL), intent (in ) :: x(:),y(size (x))
178+ real (FP_REAL) :: f(size (x))
179+ integer (FP_FLAG), optional , intent (out ) :: ierr
180+
181+ integer (FP_FLAG) :: ier
182+ integer (FP_SIZE) :: min_lwrk
183+ real (FP_REAL), allocatable :: wrk(:)
184+
185+ ! bispeu workspace: lwrk >= kx+ky+2
186+ min_lwrk = sum (this% order) + 2
187+ allocate (wrk(min_lwrk))
188+
189+ call bispeu(tx= this% t(:,1 ),nx= this% knots(1 ), &
190+ ty= this% t(:,2 ),ny= this% knots(2 ), &
191+ c= this% c, &
192+ kx= this% order(1 ),ky= this% order(2 ), &
193+ x= x,y= y,z= f,m= size (x), &
194+ wrk= wrk,lwrk= min_lwrk, &
195+ ier= ier)
196+
197+ call fitpack_error_handling(ier,ierr,' evaluate bivariate surface' )
198+
199+ end function surface_eval_many
200+
201+ ! > Evaluate surface at a single (x,y) point
202+ real (FP_REAL) function surface_eval_one(this,x,y,ierr) result(f)
203+ class(fitpack_surface), intent (inout ) :: this
204+ real (FP_REAL), intent (in ) :: x,y
205+ integer (FP_FLAG), optional , intent (out ) :: ierr
206+
207+ real (FP_REAL) :: z1(1 )
208+
209+ z1 = surface_eval_many(this,[x],[y],ierr)
210+ f = z1(1 )
211+
212+ end function surface_eval_one
213+
214+ ! > Evaluate surface on a grid domain using bispev
215+ function surface_eval_gridded (this ,x ,y ,ierr ) result(f)
216+ class(fitpack_surface), intent (inout ) :: this
217+ real (FP_REAL), intent (in ) :: x(:),y(:)
218+ real (FP_REAL) :: f(size (y),size (x))
219+ integer (FP_FLAG), optional , intent (out ) :: ierr
220+
221+ integer (FP_FLAG) :: ier
222+ integer (FP_SIZE) :: min_lwrk,min_kwrk,m(2 )
223+ real (FP_REAL), allocatable :: min_wrk(:)
224+ integer (FP_SIZE), allocatable :: min_iwrk(:)
225+
226+ m = [size (x),size (y)]
227+
228+ ! Assert real working storage
229+ min_lwrk = sum (m* (this% order+1 )) + product (this% knots- this% order-1 )
230+ if (min_lwrk> this% lwrk1) then
231+ allocate (min_wrk(min_lwrk),source= 0.0_FP_REAL )
232+ call move_alloc(from= min_wrk,to = this% wrk1)
233+ this% lwrk1 = min_lwrk
234+ end if
235+
236+ ! Assert integer working storage
237+ min_kwrk = sum (m)
238+ if (min_kwrk> this% liwrk) then
239+ allocate (min_iwrk(min_kwrk),source= 0_FP_SIZE )
240+ call move_alloc(from= min_iwrk,to = this% iwrk)
241+ this% liwrk = min_kwrk
242+ end if
243+
244+ call bispev(tx= this% t(:,1 ),nx= this% knots(1 ), &
245+ ty= this% t(:,2 ),ny= this% knots(2 ), &
246+ c= this% c, &
247+ kx= this% order(1 ),ky= this% order(2 ), &
248+ x= x,mx= size (x), &
249+ y= y,my= size (y), &
250+ z= f, &
251+ wrk= this% wrk1,lwrk= this% lwrk1, &
252+ iwrk= this% iwrk,kwrk= this% liwrk, &
253+ ier= ier)
254+
255+ call fitpack_error_handling(ier,ierr,' evaluate gridded surface' )
256+
257+ end function surface_eval_gridded
258+
147259 elemental subroutine surf_destroy (this )
148260 class(fitpack_surface), intent (inout ) :: this
149261 integer :: ierr
@@ -424,4 +536,10 @@ real(FP_REAL) function surface_derivatives_one(this,x,y,dx,dy,ierr) result(f)
424536 end function surface_derivatives_one
425537
426538
539+ ! Return fitting MSE
540+ elemental real (FP_REAL) function surface_error(this)
541+ class(fitpack_surface), intent (in ) :: this
542+ surface_error = this% fp
543+ end function surface_error
544+
427545end module fitpack_surfaces
0 commit comments