|
1 |
| -export joLinInterp1D |
2 |
| - |
3 | 1 | # 1D Linear interpolation operator
|
4 |
| -# |
5 |
| -# Parameters: |
6 |
| -# xin - input grid |
7 |
| -# xout - output grid |
8 |
| -# T - vector data type |
9 |
| -# |
10 |
| -# Note: the interval [minimum(xout),maximum(xout)] must be contained |
11 |
| -# in the interval [minimum(xin),maximum(xin)] |
12 |
| -function joLinInterp1D(xin,xout,T) |
13 |
| - xin = sort(xin) |
14 |
| - xout = sort(xout) |
15 |
| - nin = length(xin) |
16 |
| - nout = length(xout) |
17 |
| - I = nearestnbr1d(xin,xout) |
18 |
| - I[xin[I].>xout] = I[xin[I].>xout]-1 |
19 |
| - xin = [xin;1e12] |
20 |
| - b = (xout-xin[I])./(xin[I+1]-xin[I]) |
21 |
| - a = (xin[I+1]-xout)./(xin[I+1]-xin[I]) |
22 |
| - A = sparse(1:nout,I,a,nout,nin) + sparse(1:nout,min.(I+1,nin),b,nout,nin) |
23 |
| - return joMatrix(A,DDT=T,name="joLinInterp1D") |
24 |
| -end |
25 | 2 |
|
26 |
| -# 1D nearest neighbourhood operator, returns a vector I of size length(xout) |
27 |
| -# such that for every i in 1:length(xout), I[i] is the index such that |
28 |
| -# |xout[i] - xin[I[i]]| <= |xout[i]-xin[j]| for all indices j |
29 |
| -# |
30 |
| -# Parameters: |
31 |
| -# xin - input grid |
32 |
| -# xout - output grid |
33 |
| -# |
34 |
| -# Note: the interval [minimum(xout),maximum(xout)] must be contained |
35 |
| -# in the interval [minimum(xin),maximum(xin)] |
| 3 | +# helper module |
| 4 | +module joLinInterp1D_etc |
| 5 | +""" |
| 6 | +1D nearest neighbourhood operator, returns a vector I of size length(xout) |
| 7 | +such that for every i in 1:length(xout), I[i] is the index such that |
| 8 | + |xout[i] - xin[I[i]]| <= |xout[i]-xin[j]| for all indices j |
| 9 | +
|
| 10 | +Parameters: |
| 11 | + xin - input grid |
| 12 | + xout - output grid |
| 13 | +
|
| 14 | +Note: the interval [minimum(xout),maximum(xout)] must be contained |
| 15 | +in the interval [minimum(xin),maximum(xin)] |
| 16 | +
|
| 17 | +""" |
36 | 18 | function nearestnbr1d(xin,xout)
|
37 | 19 | ((minimum(xin) <= minimum(xout)) && (maximum(xin) >= maximum(xout))) || throw(ArgumentError("interval defined by xout must be contained within xin"))
|
38 | 20 | Jin = sortperm(xin)
|
@@ -63,3 +45,42 @@ function nearestnbr1d(xin,xout)
|
63 | 45 | I[Jout_inv] = I
|
64 | 46 | return I
|
65 | 47 | end
|
| 48 | +end |
| 49 | +using .joLinInterp1D_etc |
| 50 | + |
| 51 | +export joLinInterp1D |
| 52 | +""" |
| 53 | +
|
| 54 | + julia> joLinInterp1D(xin,xout,T) |
| 55 | +
|
| 56 | +1D Linear interpolation operator |
| 57 | +
|
| 58 | +# Arguments |
| 59 | +
|
| 60 | +- `xin` - input grid |
| 61 | +- `xout` - output grid |
| 62 | +- `T` - vector data type |
| 63 | +
|
| 64 | +# Signature |
| 65 | +
|
| 66 | + joLinInterp1D(xin,xout,T) |
| 67 | +
|
| 68 | +# Notes |
| 69 | +
|
| 70 | +- The interval [minimum(xout),maximum(xout)] must be contained in the interval [minimum(xin),maximum(xin)] |
| 71 | +
|
| 72 | +""" |
| 73 | +function joLinInterp1D(xin,xout,T) |
| 74 | + xin = sort(xin) |
| 75 | + xout = sort(xout) |
| 76 | + nin = length(xin) |
| 77 | + nout = length(xout) |
| 78 | + I = joLinInterp1D_etc.nearestnbr1d(xin,xout) |
| 79 | + I[xin[I].>xout] = I[xin[I].>xout].-1 |
| 80 | + xin = [xin;1e12] |
| 81 | + b = (xout-xin[I])./(xin[I.+1]-xin[I]) |
| 82 | + a = (xin[I.+1]-xout)./(xin[I.+1]-xin[I]) |
| 83 | + A = sparse(1:nout,I,a,nout,nin) + sparse(1:nout,min.(I.+1,nin),b,nout,nin) |
| 84 | + return joMatrix(A,DDT=T,name="joLinInterp1D") |
| 85 | +end |
| 86 | + |
0 commit comments