Exercise 4: Interpolation 
In this Exercise you first get a chance to play with various types of onedimensional (1D) and twodimensional (2D) splines, to get a feeling for how the splines react to motions of the knots. The interactive examples are implemented with a simple form of IDL mouse tracking, and you are invited to inspect and modify the IDL code; there are many circumstances when a quickanddirty IDL experiment or animation can be very informative.
The Problem this week consists of a (bugfree this time) set of Fortran90 procedures capable of iterative table improvement; i.e., seeking out the minimum resolution needed for a 2D table lookup to given precision. All that is needed is for you to add a 2D spline routine, patterned after the corresponding IDL routine.

Please note that the intention with bonus points is that they should be used to buy hints if you get stuck. Of course it's a fun sport to try to avoid using any, and thus maximize the score, but don't forget that the ultimate goal is to learn as much as possible. 
Subsections  approx time 

Getting the exercise files  1 min 
Play with 1D splines  60 min 
Play with 2D splines  30 min 
Write your own Fortran90 2D spline procedure  30 min 
Debug and test the program  60 min 
Home Work  2 hours 
Getting the exercise files 
[about 1 minute] 
To update the exercise files (If we have made any changes since your first CVS download), do
cd ~/ComputerPhysics cvs update d
In case of problems, see the CVS update help page.
Play with 1D splines 
[about 60 minutes] 
To start the splines.pro procedure, which demonstrates 1D splines, change directory to ComputerPhysics/4_Interpolation, and start IDLDE. To allow IDLDE to find all the needed routines that are located in the exercise directory you have to set one and possible both of the following system settings. In the new idlde version 7.0.6 you only need to set Change directory on file open under the IDL > Editor prefeerences. In 7.0.0 this does not work and you instead need to include the present work directory in the idl search path. This is done by on the IDL > Paths preference page where you inset the present exercise directory and tick it active before clicking on Apply and OK.
Then type
IDL> splines
in the command window panel.
The "splines" procedure opens an interactive plot window, with a plot of a function that is initially constant. The function is a spline that passes through a small number of points that may be moved interactively.
The type of spline is either one that is exact for 3rd order functions, or a "natural spline" (with vanishing 2nd derivative at the end points), or a local spline, with derivatives computed from the nearest neighbours, or a monotonized version of the local spline.

Play around to your hearts content, but then spend a little time playing with the IDL source code as well. First, a trivial change:
NOTE: The psym= is used in two places; first to make the initial plot, and then to update the plot (inside the loop)  you need to change in both places.

Now, let's take a look at a characteristic property of 3rd order, nonlocal splines; the influence of a single function value on the derivatives at all other points. Since the spline problem is linear (the derivatives are the solutions of a linear equation system), it is enough to play with a single point, and look at the effect on all of the derivatives. Already from the plot one can see, that setting a single function value different from the rest causes a "ringing" in the derivative values; the curve is forced to wiggle up and down to pass through all the points with constant values. The derivatives have signs that alternate, and absolute values that decay away from the perturbed point. How fast does the derivative decay? In order to find out:
repeat
loop (towards the very end of the file  for example on the cursor
command), where the derivative
is available in the variable d.
How much does the derivative decay between two points? Credits: 2/2
OK
Looking through the Section 4 lecture notes, the relative size is Credits: 2/2
OK
So, is the decay rate equal to Credits: 2/2
OK
How many intervals does it take for the ringing to drop with a factor of about 10^6 Credits: 2/2
OK
Play with 2D splines 
[about 30 minutes] 
In addition to the splines.pro procedure, there is also a similar one for 2D splines. Start this by
IDL> splines2d
and set values by clicking in the left of the two windows. This procedure has more to draw, so click, rather than drag.
Study what happens nearby a "diagonal mountain range":
IDL> device, decompose=0To find out how set this at startup, look at the Absalon IDLDE startup tutorial.

Your own Fortran90 2D spline procedure 
[about 30 minutes] 
Write a Fortran function that implements 2d interpolation (cf. the lecture notes), either by translating the IDL routine spl2d.pro, or from scratch. This should take less than 30 minutes.

You may either choose to use the idl2f90 filter, with
idl2f90 < spl2d.pro > spl2d.f90and then clean up the result, or else (probably more efficient as a learning exercise) write it from scratch (it is still OK to cut & paste the complicated expressions from spl2d.pro).
Overall structure of the function:
Module spl2df contains function spl2d (....,nx,ny) ! LOOK AT THE CALL (nx and ny are there!) implicit none ! forces everything to be defined explicitly integer:: ... ! integer variables real:: ... ! real variables, some with dimensions ! compute indices to the four nearest corners from x,y ... [ Here things depend on if you dimension to (nx,ny) or (0:nx1,0:ny1) ! put limits on the integer addresses ... [ use Fortrans min() and max() functions ] ! compute fractions px, py ... [ make sure that these allow px larger than 1 when x is outside ] ! compute g and h coeffs in all four corners ... [ Use cut & paste to take this from the IDL function ] ! compute the interpolated value ... [ Cut & paste this too from spl2d.pro ] ! return the value ... [ A couple of changes are needed to do it the Fortran way ] end function spld2d end module
Need a hint or four?




So, based on the careful preparations you did before you came (:), you pretty soon have something you think is a correct function. Well, make sure there are at least no syntax errors, by compiling the function:
localhost:home> cd ComputerPhysics/4_Interpolation localhost:4_Interpolation> gfortran c spl2d.f90
Debug and test the program 
[about 60 minutes] 
So, how do you know that the program works (apart from a false feeling of confidence :)? The hint is free of charge: Testing!
There is a small test program, Test.f90, that tests your interpolation routine on simple polynomials, where the result should be exact. In "real life", you should write such test programs your self; they are typically a smaller effort to code than the program piece they can test and thus, with a small extra investment, you obtain a very valuable result; "true" ^{1)} confidence that the code works.
^{1)} Well, there are always catches; all a test really proves is that the code works in the test situation. Whether or not one can be confident that it also works in the application depends on how similar the situation is (i.e., on how good  and bug free  the test program is :).
Anyway, in order to compile the subroutine and the test program, and run the test, here is what you do:
localhost:4_Interpolate> gfortran Test.f90 spl2d.f90
localhost:4_Interpolate> ./a.out
The program tries successively 1st, 2nd, and 3rd order functions, prints an error estimate followed by OK or NOK. The test covers both the case where all points lie inside (interpolations proper), and the case where some points lie outside (extrapolations). The latter case is there to check if the bounds on the indices are correct.
If it still does not work, you need to use the DDD debugger, but even if it does, you should try it out in this simple situation:
gfortran
compiler and a g
option:
localhost:4_Interpolate> gfortran g spl2d.f90 Test.f90
localhost:4_Interpolate> ddd a.out & < free hint here!!
The testing phase should leave you with a spl2d that you may be confident works, also in the full program.
To compile and run the iterative table improvement program Lookup.x:
localhost:4_Interpolate> make clean localhost:4_Interpolate> make
The first command removes the output from the gfortran
compilations, to
allow the second "make" to work also at fys.ku.dk, where only the ifort
compiler is available.
You will learn more about the make
command later, for
now it
is enough to know that, using a description of the relation between
pieces
of a program, stored in a file named Makefile, it may be
used to
simplify the compilation of programs with more than just one or two
subroutines.
Here, it has been setup to produce the executable file Lookup.x
and then
run it once.

localhost:4_Interpolate> ./Lookup.x
After you are done with the exercise, clean up the directory, to save space. Remove *.o, *.x, idl.ps, and a.out files. A convenient way of doing this is built into the Makefile; the target clean does the cleanup, so type
localhost:4_Interpolate> make clean
Upload the spl2d.f90 file below.

Your home work 
[about 2 hours] 
If you don't feel up to date with Fortran, then spend time reading in your Fortran book, about declarations, subroutines and functions, and about arrays and array indices. You should also follow up on the last lecture, by browsing the Chapter about interpolation in Numerical Recipes. Note, however, that the spline discussion there is unnecessarily complicatedconsider it cursory material. As a preparation for the next Lecture, it is also a good idea to browse the Chapter in Numerical Recipes about numerical integration (quadrature).