;--------------------------------------------------------------------
function test_idl_interpos_all_cases, case_to_test, tension_in, option_in
; make interposgfortran_idl (to create the .so library, in ../interpos_libs folder)

  libname='/afs/ipp-garching.mpg.de/home/o/osauter/interpos_develop/interpos_libs/gfortran_idl/libinterposgfortran_idl.so'
; on toks0 machines:
  ;libname='/afs/ipp-garching.mpg.de/home/o/osauter/bin/toks0/libinterposgfortran_idl.so'
; on sxaug machines (still ELF error)
  ;libname='/afs/ipp-garching.mpg.de/home/o/osauter/bin/sxaug/libinterposgfortran_idl.so'
   libname='/home/sauter/bin/libinterposgfortran_idl.so' ; on lac3 at spc
   libname='../interpos_libs/gfortran_cwrapper/libinterposgfortran_cwrapper.so' ; on lac3 at spc
; basic case of a profile
  ioptout = long(1)
  nout = long(101)
; for profile cases, basic data points
  xin=double([0.1, 0.2, 0.3, 0.5, 0.7, 1.0])
  yin=double([1.3, 1, 0.6, 0.3, 0.2, 0.05])
  nin = long(n_elements(xin))
  nbc = long([1 ,0])
  ybc = double([0., 0.])
  info = long(-1)
  xout = 0 + DINDGEN(nout)/(nout-1)*(xin[nin-1]-0)
; for periodic cases, basic data
  pi = double(3.141592653589793)
  ninper = long(61)
  t=FINDGEN(fix(ninper))/(ninper-1)*2*pi
  xinper=double(t)
  yinper = double(sin(0.3+xinper)+0.9*randomu(seed,fix(ninper)))
  nbcper = long(-1)
  ybcper = double(2*pi)
  xoutper = double((DINDGEN(fix(nout))/(nout-1)*4-1)*pi)
; need to create memory space for outputs
  yout=DblArr(fix(nout))
  youtp=yout
  youtpp=yout
  youtint=yout
  if (n_elements(tension_in) EQ 0) THEN BEGIN
     tension=double(-0.1)
  ENDIF ELSE BEGIN
     tension=double(tension_in)
  ENDELSE
  if (n_elements(option_in) EQ 0) THEN BEGIN
     option = long(63)
  ENDIF ELSE BEGIN
     option=long(option_in)
  ENDELSE
  
;; call external fitting routine
  s=999
  CASE case_to_test OF
     1: BEGIN
        print, 'case 1, profil with only yout and basic bnd conditions. Try: aa=test_idl_interpos_all_cases(1,-0.01)'
; standard call but extrapolating towards 0
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin, yin, nin, yout, tension, xout, nout)
        window, 0
        plot,xin,yin,yrange=[0, 2],psym=4
        oplot,xout,yout,linestyle=1
        xin_withaxis = DblArr(n_elements(xin)+1)
        xin_withaxis[1:nin] = xin
        xin_withaxis[0] = 0
        yin_withaxis = DblArr(n_elements(yin)+1)
        yin_withaxis[1:nin] = yin
        yin_withaxis[0] = yin[0] ; do not know how else to determine value on axis at this stage
; this way point at left is at rho=0 and we can impose 1st derivative=0
        nbc = long([1 ,0])
        ybc = double([0., 0.])
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin_withaxis, yin_withaxis, nin+1, yout, tension, xout, nout, nbc, ybc)
        oplot,xout,yout,COLOR='FFF000'x,linestyle=2
; add error bar to leave first point free
        sigma = DblArr(nin+1)
        sigma[*]=double(1.)
        sigma[0] = double(100.)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin_withaxis, yin_withaxis, nin+1, yout, tension, xout, nout, nbc, ybc, sigma)
        oplot,xout,yout,COLOR=160
        return,yout
        tension = double(0.)    ; compare with standard cublic spline (thus 2nd der=0 and no extra central point)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin, yin, nin, yout, tension, xout, nout)
        oplot,xout,yout,linestyle=3
        print,' diamonds=data,blackdotted=tens_extrap,balckdashdot=tension=0,bluedash=1stder=0no_errorbar, red=1stder=0_sigma'
     END
     2: BEGIN
        print, 'as case 1 but tests derivatives and other inputs, only from best fit with sigma and 1st der=0 in center'
        xin_withaxis = DblArr(nin+1)
        xin_withaxis[1:nin] = xin[*]
        xin_withaxis[0] = 0
        yin_withaxis = DblArr(n_elements(yin)+1)
        yin_withaxis[1:nin] = yin
        yin_withaxis[0] = yin[0] ; do not know how else to determine value on axis at this stage
        sigma = DblArr(nin+1)
        sigma[*]=double(1.)
        sigma[0] = double(100.)
        nbc = long([1 ,0])
        ybc = double([0., 0.])
        ioptout = long(2)       ;
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin_withaxis, yin_withaxis, nin+1, yout, tension, $
                          xout, nout, nbc, ybc,sigma,youtp)
        window, 0
        plot,xin,yin,psym=4
        oplot,xout,yout,COLOR=160
        window,1
        plot,xout,-youtp/(abs(yout)+3e-2),COLOR=160,title='scalelength = -1st derivative/function'
        print,' in wshow,1: -1st der/function derivative (scalelength)'
        ioptout = 3L            ;
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin_withaxis, yin_withaxis, nin+1, yout, tension, $
                          xout, nout, nbc, ybc, sigma,youtp,youtpp)
        wset,1
        oplot,xout,youtp,linestyle=2

        window,2
        plot,xout,youtpp,title='2nd derivative'
        print,' in wshow,2: 2nd derivative'
        ioptout = 4L            ;
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin_withaxis, yin_withaxis, nin+1, yout, tension, $
                          xout, nout, nbc, ybc, sigma,youtp,youtpp,youtint)
        wset,2
        oplot,xout,youtpp,linestyle=2

        window,3
        plot,xout,youtint,title='integral'
        print,' in wshow,3: integral'
     END
     3: BEGIN
        print, 'as case 1 but force value at axis and edge'
        nbc = long([2 ,2])
        ybc = double([1.6, 0.02])
        ioptout = long(1)       ;
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin, yin, nin, yout, tension, xout, nout, nbc, ybc)
        window, 0
        plot,xout,yout,COLOR=160 ; this is very wrong, because value at xin[0]=1.6, but xin[0]=0.1 not 0, so do it at 0 by adding a point
        oplot,xin,yin,psym=4
        xin_withaxis = DblArr(nin+1)
        xin_withaxis[1:nin] = xin[*]
        xin_withaxis[0] = 0
        yin_withaxis = DblArr(n_elements(yin)+1)
        yin_withaxis[1:nin] = yin
        yin_withaxis[0] = yin[0] ; exact value not important here since set by ybc
        s = call_external(libname,'interpos_wrapper_c',ioptout, xin_withaxis, yin_withaxis, nin+1, yout, tension, xout, nout, nbc, ybc)
        oplot,xout,yout,COLOR='FFF000'x
        print,' value at left set to 1.6 and at right to 0.02. Red, left point is 1st data point at rho=0.1'
     END
     4: BEGIN
        print, 'fit and extrapolations'
        ioptout = long(1)
; minimal call for fit on same xout=xin
        yout=DblArr(fix(ninper))
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout)
        window, 0
        plot,xinper,yinper,psym=4,xrange=[0.,2*pi],xstyle=1,title='interpolation of sin+random'
        oplot,xinper,yout,COLOR=160,thick=2
; minimal call with given tension
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension)
        oplot,xinper,yout,COLOR='FFF000'x,thick=2, linestyle=3
; interpolate and extrapolate on xout
        yout=DblArr(nout)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout)

        window, 1
        plot,xoutper,yout,COLOR='FFF000'x,thick=2, linestyle=0,title='extrapolation of sin(0.3+x)+random',yrange=[-3.,3.],ystyle=1
        oplot,xinper,yinper,psym=4
; add option for extrapolation, need to fill in default values for other parameters before option
        nbc=long([0, 0])
        ybc=double([0., 0.])
        sigma = DblArr(ninper)
        sigma[*] = double(1.)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout, $
                          nbc,ybc,sigma,option)
        oplot,xoutper,yout,COLOR=160,thick=2
        print,'change 3rd input argument ''option'' to change extrapolation function, see http://spc.epfl.ch/interpos (scroll down)'
        print,'last digit of option, 1, 2, or 3 means linear, quad or cubic interpolation, 1st digit and sign yields extrapolation'
        return,yout
     END
     5: BEGIN
        print, 'same as 4 but use periodic boundary conditions'
        ioptout = long(1)
; minimal call for fit on same xout=xin
        yout=DblArr(fix(ninper))
        nbc[0] = long(-1)
        ybc[0] = double(2*pi)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout,tension,xinper,ninper,nbc,ybc)
        window, 0
        plot,xinper,yinper,psym=4,xrange=[0.,2*pi],xstyle=1,title='interpolation of sin+random with periodic b.c'
        oplot,xinper,yout,COLOR=160,thick=2
; obtain the same result if nbc, ybc have length(1), since only 1st value used if nbc(1)=-1
        nbc = long(-1)
        ybc = double(2*pi)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout,tension,xinper,ninper,nbc,ybc)
        oplot,xinper,yout,COLOR='FFF000'x,thick=2,linestyle=3

; interpolate and extrapolate on xout (note extrapolation option is not used since "extrapolates" with periodic
; conditions, but last digit used for interpolation. Periodic conditions only for cubic spline
        yout=DblArr(nout)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout,nbc,ybc)

        window, 1
        plot,xoutper,yout,COLOR='FFF000'x,thick=2, linestyle=0,title='extrapolation of sin(0.3+x)+0.3random',yrange=[-2.,2.],ystyle=1
        oplot,xinper,yinper,psym=4
; add option for interpolation, need to fill in default values for other parameters before option
        sigma = DblArr(ninper)
        sigma[*] = double(1.)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout, $
                          nbc,ybc,sigma,option)
        oplot,xoutper,yout,COLOR=160,thick=2
        print,'interpolation and extrapolation performed with periodicity 2pi condition, thus option not used with cubic spline'
        print,'last digit of option, 1, 2, or 3 means linear, quad or cubic interpolation, 1st digit and sign yields extrapolation'
        return,yout
     END
     6: BEGIN
        print, 'compute interpolation/extrapolation with periodic boundary conditions and derivatives'
        yinper = double(sin(xinper)+0.02*randomu(seed,fix(ninper)))
        ioptout = long(1)
        yout=DblArr(nout)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout,nbcper,ybcper)
        window, 0
        plot,xoutper,yout,COLOR=160,thick=2, linestyle=0,title='extrapolation of sin(x)+0.02random',yrange=[-2.,2.],ystyle=1
        oplot,xinper,yinper,psym=4
        youtp=DblArr(nout)
        youtpp=DblArr(nout)
        youtint=DblArr(nout)
        ioptout = long(2)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout, $
                          nbcper,ybcper,youtp)
        ioptout = long(3)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout, $
                          nbcper,ybcper,youtp,youtpp)
        ioptout = long(4)
        s = call_external(libname,'interpos_wrapper_c',ioptout, xinper, yinper, ninper, yout, tension, xoutper, nout, $
                          nbcper,ybcper,youtp,youtpp,youtint)
        window, 1
        plot,xoutper,youtp,COLOR='FFF000'x,thick=2,linestyle=0,title='1st derivative, should be close to cos(x)'
        window, 2
        plot,xoutper,youtpp,COLOR='FFF000'x,thick=2,linestyle=1,title='2nd derivative, should be close to -sin(x)'
        window, 3
        plot,xoutper,youtint,COLOR='FFF000'x,thick=2,linestyle=2,title='integral from x=0, should be close to 1-cos(x)'
     END
  END
  return, info
end