I took @urbanjost 's “inputfile”,
prefixed it with ncolumn and nrowdata ,
appended C0,T0 values list.
There is considerable non-linearity in the surface generated, especially near cm=0.1, tmu=15. which shows that linear interpolation is poor for the saddles (see below)
However, I then included a bi-linerar interpolation estimate to the program, based on the 4 surounding table values and provided an estimate.
Estimates near the identified saddle show significant variation are probably not adequate, but the program as presented does achieve values, hopefully better than the previous nearest table samples.
The file ‘muHA.dat’ is
5 7
0 .1 .2 .3 .4 .5
5 20526.3716 14572.9521 27718.159 22593.21442 27145.10622
10 32401.4931 5641.11114 27078.18790 19103.3273 5604.23500
15 855.775 20600.1222 25916.23962 25393.20786 23221.16154
20 29004.23590 20097.9097 12671.16411 17876.26376 30291.635
25 8187.20373 22478.21666 1628.3514 19902.32148 8906.4292
30 16514.10185 32576.30782 23101.1582 21897.16568 8106.30035
35 6549.6536 26006.25288 19449.24701 14215.14241 21347.9012
.2 25
.25 25
.25 27
.12 14.
.12 15.
.12 16.
-1 -1
The program, including a bilinear interpolation is now:
implicit none
integer, parameter :: wp = selected_real_kind (10)
real(wp), external :: get_linear_value
real(wp), allocatable :: cm(:), Tmu(:), amur(:,:)
real(wp) :: T0, C0, val_amur
integer :: ncm, ntmu, kc, kt, i, j, iostat
ncm = 11
ntmu = 269
T0 = 300.0_wp
c0 = 0.4_wp
kc = 1 ! closest cm(i) for c0
kt = 1 ! closest Tmu(jt) for T0
open(10,file='muHA.dat')
read (10,*) ncm, ntmu
write (*,*) ncm, ' columns of cm values'
write (*,*) ntmu,' rows of Tmu values'
allocate ( cm(ncm), Tmu(ntmu), amur(ncm,ntmu) )
read (10,*, iostat = iostat) Tmu(1), (cm(i),i=1,ncm)
if ( iostat /= 0 ) then
write (*,*) ' problem reading cm from muHA.dat'
stop
end if
do j = 1,ntmu
read (10,*, iostat = iostat) Tmu(j), (amur(i,j),i=1,ncm)
if ( iostat /= 0 ) then
write (*,*) ' problem reading amur from muHA.dat at record',j
exit
end if
end do
if ( j <= ntmu ) then
print*, 'unable to read table'
stop
else
write (*,*) ntmu,' rows of amur values read'
end if
do
read (10,*,iostat=iostat ) C0,T0
if ( iostat /= 0 ) exit
write (*,*) ' '
write (*,*) 'Testing for C0=',C0,' T0=',T0
do i=1,ncm
if ( abs(cm(i)-c0) < abs (cm(kc)-c0) ) kc = i
end do
do j = 1,ntmu
if ( abs(Tmu(j)-T0) < abs(Tmu(kt)-T0) ) kt = j
end do
write (*,*) 'Closest value for'
write (*,*) ' C0 =',c0,' is', cm(kc)
write (*,*) ' T0 =',t0,' is', Tmu(kt)
write (*,*) ' amur in table =',amur(kc,kt)
val_amur = get_linear_value ( C0, T0, cm, Tmu, amur, ncm, ntmu )
write (*,*) ' linear interpolation amur =',val_amur
end do
close (10)
end
function get_linear_value ( C0, T0, cm, Tmu, amur, ncm, ntm )
! does bi-linear interpolation
! assumes cm and Tmu are ordered; reports error if out of order
! if C0 or T0 are out of range; limits to edge of table
implicit none
integer, parameter :: wp = selected_real_kind (10)
real(wp) :: get_linear_value
integer :: ncm, ntm
real(wp) :: cm(ncm), Tmu(ntm), amur(ncm, ntm)
real(wp) :: T0, C0
!
integer :: kc, kt
real(wp) :: fc, ft
call get_linear_factor ( fc, kc, C0, cm, ncm, 'cm' )
if ( kc < 1 ) goto 901
call get_linear_factor ( ft, kt, T0, Tmu, ntm, 'Tmu' )
if ( kt < 1 ) goto 901
write (*,*) 'interpolation factors :',kc,kt, fc,ft
!
get_linear_value = amur(kc, kt) * (1-fc)*(1-ft) &
+ amur(kc+1,kt) * fc *(1-ft) &
+ amur(kc, kt+1) * (1-fc)* ft &
+ amur(kc+1,kt+1) * fc * ft
return
901 get_linear_value = -999.
end function get_linear_value
subroutine get_linear_factor ( fc, kc, C0, cm, ncm, vn )
implicit none
integer, parameter :: wp = selected_real_kind (10)
integer :: kc, ncm
real(wp) :: fc, C0, cm(ncm), dx
character :: vn*(*)
do kc = 1,ncm-1
dx = cm(kc+1)-cm(kc)
if ( dx == 0 ) cycle ! step in values so no inetrpolation
fc = (C0-cm(kc)) / dx
if ( dx < 0 ) goto 901 ! CM out of order
if ( fc >= 0 .and. fc <= 1 ) exit
end do
if ( kc < ncm ) return
if ( C0 < cm(1) ) then
kc = 1
fc = 0
return
else if ( C0 > cm(ncm) ) then
kc = ncm-1
fc = 1
return
else
goto 902
end if
901 write (*,*) ' ERROR ',vn,' list is out of order'
kc = -1
return
902 write (*,*) ' ERROR ',vn,' target ',C0,' is out of range',cm(1),cm(ncm)
kc = -1
return
end subroutine get_linear_factor
This bilinear interpolation is fairly straightfoward; identifying the 4-point enclosing cell and factors.
The next step could be to introduce 9-point quadratic or 16-point cubic interpolation.
The original table posted also had non uniform column and row spacing for cm and Tmu, which is another issue which quadratic or cubic interpolation might ignore.
At this stage, I think it may be preferable to better understand the field behaviour or increase the table resolution.
