Linear interpolating between corresponding matrix coefficients..... 1

[aneedfortheory]

Looking for help from Fortran programmers.....
I have a problem with linear interpolation. I have constructed a subroutine that reads from an external array(s) that is called rates(tk,u,l) where tk is the temperature (that helps distinguish between the various arrays, u for the upper level (i'm dealing with transitions) and l for lower. For a temperature not mentioned in the list of temperatures the code is instructed to interpolate the tabled values to construct an array for that specific tyemperature. Therefore the routine calls each matrix coefficient in turn. The problem i have is that what if the "rate" or matrix coefficient corresponding to the higher interpolating value or temperature is below that of the coefficient corresponding to the lower value together with rates at which the situation is switched. I have the coding written below and would be grateful if someone was to follow through it and specify where i should make changes....Thanks in advance!!

Regards,
Robert Loughnane.

PS: the variable 'it' corresponds to a coefficient that cycles through the 8 temperatures in turn but only the matrix values corresponding to any two of these temperatures is used at any one time.....


subroutine colrateion(tk)

c A linear interpolation of the ion rates.
c Using the tables of temperatures and rates in
c common IONRATES, interpolates the rates for
c the temperature TK passed as an argument.
c The interpolated rates are stored in the common
c COLL.

c CR Interpolated collision rates cm^3/sec at
c temperature TK
c NSTATE number of levels
c TEMPS temperatures in K
c RATE collision rates at temperatures in TEMPS
c IT array index

real * 8 cr
integer u,l
parameter (nstate=22)
common /coll/ cr(nstate,nstate)
common /ionrates/ temps(8),rate(8,22,22)

it = tk/10. + 1
it = max(1,min(it,7))

do 2 u = 1,nstate
do 2 l = 1,nstate

c This is a simple linear interpolation.
if ((rate(it+1,u,l) - rate(it,u,l)) .lt. 0.) then
cr(u,l) = dble (rate(it,u,l)
& + (tk - temps(it))/(temps(it+1) - temps(it))
& * (rate(it,u,l) - rate(it+1,u,l)) )
else
cr(u,l) = dble (rate(it,u,l)
& + (tk - temps(it))/(temps(it+1) - temps(it))
& * (rate(it+1,u,l) - rate(it,u,l)) )
endif

write(6,*) tk,rate(it,u,l),rate(it+1,u,l),cr(u,l),u,l

2 continue

return
end

 

[xwb]

1) If it was outside the range then you'd probably have to extrapolate the values. Don't really know much about extrapolation, also it depends on how far out they are and whether it is linear or not. THere are lots of algorithms for extrapolation using difference methods.

2) You could reduce the code by using abs. i.e.

              cr(u,l) = dble (rate(it,u,l)
     &        + (tk - temps(it))/(temps(it+1) - temps(it))
     &        * abs(rate(it+1,u,l) - rate(it,u,l)) )

Then you wouldn't need the if statement.

3) Are you Dr Robert from Aston who used to work with David Scrimshire?

 

[aneedfortheory]

xwb,
Thanks for the prompt answer, i have tried abs as well as fabs and the way in which you constructed the above bit of code was giving a little trouble. Extrapolation will eventually be of a concern, but as of yet, only interpolation matters much. I'm not exactly that Dr. Robert, i'm his 23 year old eldest son...he studied Biomedical Engineering, i'm a Molecular Astrophysicist!

Robert.

PS: By the way, is there any other more shrewd, not too wayward, that can interpolate bewteen coefficients by taking into account the upness and downness of the alrenate values?

 

[xwb]

One way I can think of is a linear regression curve fit and to find the point on the curve but for the size of matrix you are dealing with, it could be quite expensive in time.

Another way would be to get the values between it and it+1, (it-1 and it+1) and (it and it+2) and average them.

Regards to your dad - you could ask him about extrapolation methods. His first degree was Maths.