# document.write(pagetitle)

Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5
```;+
; NAME:           ULY_SSP_INTERP
;
; PURPOSE:        Get interpolated spectrum for logT, Met and [Mg/Fe]
;
; USAGE:          model_new = uly_ssp_interp(model_grid, para)
;
; DESCRIPTION:    Interpolate spectrum from a grid of models at given
;                 log(age), [Fe/H] and [Mg/Fe]. Use ULY_2DERIV and
;                 INTERP_3D. Called, for example, from ULY_SSP_EXTR
;                 and ULY_FIT_LIN.
;
; ARGUMENTS:
;   model_grid:   Grid of model (structure)
;   [x, y, mgfe]: Log of age (Myr); metallicity (dex); [Mg/Fe]
;
; RETURN:         Interpolated spectrum (array)
;
; EXAMPLE:        Read a model grid and then extract a spectrum at
;                 given 10 Gyr and solar metallicity
;                 model = uly_ssp_interp(model_grid, [alog(10000), 0.0])
;
; AUTHOR:         Philippe Prugniel
;-
; CATEGORY:       ULY_SSP
;-----------------------------------------------------------------------------

function uly_2deriv,x,y
;
; NAME:         ULY_2DERIV
;
; PURPOSE:      Compute the 2nd derivatives of Y(L,N) at points Xi (i=0,..N)
;
; USAGE:
;               y2 = uly_2deriv(x,y)
; DESCRIPTION:
;               Given arrays X of length N and Y of size LxN containing a
;               tabulated function, i.e.
;               Yi = f(Xi), with X1 < X2 < ... < Xn,
;               this routine returns and array Y2 of size LxN which contains
;               the second derivatives in the second coordinates of the
;               interpolating function at the tabulated points Xi.
;
;               Used by the functions ULY_SSP_INTERP and DERIVE_3D
;
;
; ARGUMENTS:
;           x:  independent variable vector
;           y:  dependent variable vector
;
; RETURN:
;           y2: second derivatives at all x, of length n
;               the type of y2 is of the highest precision of x and y
;
; HISTORY:
;      D. Neill, October, 1991
;      http://www.astro.washington.edu/deutsch-bin/getpro/library43.html?SPLINT
;      adapted for ULySS,  Mina Koleva & Philippe Prugniel 20070308
;-
;------------------------------------------------------------------------------

s = size(y)
n = s[2]

if s[2] ne n_elements(x) then message, 'X and Y should have the same length'

type = size(y, /TYPE) > size(x, /TYPE)
y2 = make_array(s[1],n, /NOZERO, TYPE=type)
u = make_array(s[1],n, /NOZERO, TYPE=type)
;
; The lower boundary condition is set to be "natural"
;
y2[*,0] = 0.
u[*,0] = 0.

;
; This is the decomposition loop of the tridiagonal algorithm.  Y2 and
; U are used for temporary storage of the decomposed factors.
; We make the choice of using loops along the interpolating region
; rather than using TRISOL or SPL_INIT, because this allows to use
; array operations along the other axis. This is the best solution because
; this other axis is in principle much longer (i.e. it is the number
; of wavelength points and is >10 times larger than the length of the
; interpolated axis).
;
tsig = double((x - shift(x, -1))) / (shift(x, +1) - shift(x, -1))

status = check_math(MASK=32)  ; this function may generate underflows

for i=1,n-2 do begin
p = tsig[i] * y2[*,i-1] + 2.
y2[*,i] = (tsig[i]-1.) / p

u[*,i]=( 6. * ((y[*,i+1]-y[*,i]) / (x[i+1]-x[i]) - (y[*,i]-y[*,i-1]) \$
/ (x[i]-x[i-1])) / (x[i+1]-x[i-1]) - tsig[i]*u[*,i-1]) / p
endfor

; The upper boundary condition is set either to be "natural"
;
y2[*,n-1] = 0.
;
; This is the backsubstitution loop of the tridiagonal algorithm
;
for k=n-2,0,-1 do begin
y2[*,k] = y2[*,k] * y2[*,k+1] + u[*,k]
endfor

if status eq 0 then status = check_math(MASK=32)  ; dismiss underflows

return, y2
end

;-----------------------------------------------------------------------------
function interp_3d, x, y, TMPARR=tmparr, D2X=d2x,\$
XGRID=xgrid, YGRID=ygrid
compile_opt idl2, hidden

; NAME:         INTERP_3D (internal routine)
;
; PURPOSE:      Interpolate in a data cube
;
; USAGE:        z = interp_3d( x, y, TMPARR=tmparr, D2X=d2x,\$
;                               XGRID=xgrid, YGRID=ygrid)
; DESCRIPTION:  Interpolate in a data cube (tmparr) at (x,y) coordinates
;
; ARGUMENTS:
;      x, y:    The demanded coordinates at which the cube should be
;               intepolated
;
; KEYWORDS:
;   TMPARR :    model cube
;   D2X    :    array of X 2nd derivatives (input)
;   XGRID  :    array of X for the cube (input)
;   YGRID  :    array of Y for the cube (input)
;
; HISTORY:       <unknown> author

if not finite(x) then message, 'Input (X) is not finite'

s = size(tmparr)

if(x gt max(xgrid)) then x = max(xgrid)
if(y gt max(ygrid)) then y = max(ygrid)
if(x lt min(xgrid)) then x = min(xgrid)
if(y lt min(ygrid)) then y = min(ygrid)

if(s[0] eq 1) then return, tmparr

;    t=systime(1)

if((s[0] eq 3) and keyword_set(d2x) and (s[3] ge 4)) then begin

; determine in what region is the point[x,y]
xval = interpol(findgen(s[2]), xgrid, x) ;
yval = interpol(findgen(s[3]), ygrid, y)
xmin = floor(xval)
xmax = xmin + 1
if(xmax ge s[2]) then xmax=xmin
ymin = floor(yval)
ymax = ymin + 1
if(ymax ge s[3]) then ymax=ymin
ys1 = ymin - 1              ;limit for the spline
if(ys1 lt 0) then ys1 = ys1 + 1
ys2 = ys1 + 3    ;we need 4 points of metallicities to make spline
if(ys2 ge s[3]) then begin
ys2 = s[3] - 1
ys1 = ys2 - 3
endif

intvec = dblarr(s[1], /NOZERO) ;init. interpolated vector

;interpolate along x direction for the 4 y nodes that we need for the
;spline in y

h = (xgrid[xmax]-xgrid[xmin])
ax = (h eq 0d) ? 0 : (xgrid[xmax]-x) / h
;    ax = (xgrid[xmax]-x) / (xgrid[xmax]-xgrid[xmin])
bx = 1.0d - ax
cx = (ax^3-ax) * ((xgrid[xmax]-xgrid[xmin])^2) / 6d
dx = (bx^3-bx) * ((xgrid[xmax]-xgrid[xmin])^2) / 6d
mm = reform(ax*tmparr[*,xmin,ys1:ys1+3]+bx*tmparr[*,xmax,ys1:ys1+3]+\$
cx*d2x[*,xmin,ys1:ys1+3]+dx*d2x[*,xmax,ys1:ys1+3], s[1], 4)
d2y = uly_2deriv(ygrid[ys1:ys1+3],mm)

h = (ygrid[ymax]-ygrid[ymin])
ay = (h eq 0) ? 0 : (ygrid[ymax]-y)/h
;    ay = (ygrid[ymax]-y)/h
by = 1.0d - ay
cy = ay^3 - ay
dy = by^3 - by
intvec[*] = ay*mm[*,ymin-ys1]+by*mm[*,ymax-ys1]+\$
(cy*d2y[*,ymin-ys1]+dy*d2y[*,ymax-ys1])*(1d/6d)*(h^2)
endif else begin
if s[0] ne 3 then message, 'The model grid must be 3D'
if s[3] le 4 then message, 'The grid must have > 3 points for spline interpolation'
endelse

; print, 'TIME INTERP 3D',systime(1)-t

return, intvec
end

; ----------------------------------------------------------------------------
function uly_ssp_interp, model_grid, para
compile_opt idl2, hidden; para[0] : log age
; para[1] : metallicity
; para[2] : Mg/Fe

; Diagnostic required extrapolations
m0 = min(*model_grid.o_age, MAX=m1)
if para[0] lt m0 or para[0] gt m1 then message, /INFO, 'Requested age is out of range ... extrapolation'
m0 = min(*model_grid.o_metal, MAX=m1)
if para[1] lt m0 or para[1] gt m1 then message, /INFO, 'Requested [Fe/H]='+\$
strtrim(para[1],2)+' is out of the valid range ('+strtrim(m0,2)+\$
','+strtrim(m1,2)+')  ... extrapolation'

if size(*model_grid.data, /N_DIM) gt 3 then begin  ; grid with Mg/Fe resol
if size(*model_grid.data, /N_DIM) ne 4 \$
or (size(*model_grid.data, /DIM))[3] ne 2 then \$
message, 'Support on grid with 2 [Mg/Fe] levels)'
inpmgfe = *model_grid.o_mgfe
const = (exp(inpmgfe[1])-exp(para[2])) / (exp(inpmgfe[1])-exp(inpmgfe[0]))
return, const * \$
interp_3d(para[0], para[1], TMPARR=(*model_grid.data)[*,*,*,0], \$
D2X=(*model_grid.d2t)[*,*,*,0],\$
XGRID=*model_grid.o_age, YGRID=*model_grid.o_metal) + \$
(1. - const) * \$
interp_3d(para[0], para[1], TMPARR=(*model_grid.data)[*,*,*,1], \$
D2X=(*model_grid.d2t)[*,*,*,1],\$
XGRID=*model_grid.o_age, YGRID=*model_grid.o_metal)

endif else \$
return, interp_3d(para[0], para[1], TMPARR=*model_grid.data, \$
D2X=*model_grid.d2t,\$
XGRID=*model_grid.o_age, YGRID=*model_grid.o_metal)
end

; ---- end -------------------------------------------------------------------
```
Part of ULySS package