;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO getsize, $
   series

; this subroutine determines the array sizes 
; needed to store all the data

@compile_opt

IF N_ELEMENTS(series) EQ 0 THEN series='u'

dir00='/u/alex/'
dir0=dir00+'kepler/'+series+'nuc/'

RESOLVE_ROUTINE, $
   series+"nuc", $
   /COMPILE_FULL_FILE, $
   /NO_RECOMPILE  
CALL_PROCEDURE, $
   'get'+series+'runs', $
   mass, $
   METALLICITY=metallicity, $
   EXPLOSIO=exp

dump="#nucleo"

nmass=N_ELEMENTS(mass)
nexp=N_ELEMENTS(exp)

jmmax=0
imaxbe=INTARR(nexp)

nions=0
ions=LONARR(nexp,10000)
nions=LONARR(nexp)

FOR im=0,nmass-1 DO BEGIN
    ms=mass_string(mass[im])
    name=series+ms
    dir=dir0+name+'/'
    FOR ie=0,nexp-1 DO BEGIN
        file=dir+name+exp[ie]+dump
        loaddump,file
        getdump,'nq 53',imaxb
        getdump,'nq 2',jm
        getdump,'idxb',idxb
        jmmax>=jm
        imaxbe[ie]>=imaxb
        FOR ii=0,imaxb-1 DO BEGIN
           ij=(WHERE(ions[ie,0:nions[ie]-1>0] EQ idxb[ii]))[0]
           IF ij EQ -1 THEN BEGIN
              ions[ie,nions[ie]]=idxb[ii]
              nions[ie]+=1 
           ENDIF
        ENDFOR
    ENDFOR
ENDFOR

PRINT,'jm_max         =         ',long(jmmax)
FOR ie=0,nexp-1 DO BEGIN
    PRINT,'imaxb-'+exp[ie]+'        = ',imaxbe[ie],nions[ie]
ENDFOR
PRINT,'imaxb-max      = ',MAX(imaxbe),MAX(nions)

ions_all=REFORM(ions[0,0:nions[0]-1])
nions_all=nions[0]
FOR ie=1,nexp-1 DO BEGIN
   ions_all=[ions_all,REFORM(ions[ie,0:nions[ie]-1])]
ENDFOR
ions_all=UNIQ(ions_all,SORT(ions_all))
nions_all=N_ELEMENTS(ions_all)
PRINT,'imaxb-all-ions =         ',nions_all


END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO make_nucleo_zone_alliso_bgal, SERIES=series, EXP=exp

@compile_opt

ftag = ' [MAKE_NUCLEO_ZONE_ALLISO_BGAL] '

IF N_ELEMENTS(series) EQ 0 THEN series='u'

dir00='/u/alex/'
dir0=dir00+'kepler/'+series+'nuc/'

CASE series OF
   'z':BEGIN 
      maxzone=1105
      maxiso=2750
   END
   'u':BEGIN 
      maxzone=1222
      maxiso=1842
   END
   ELSE: BEGIN
      maxzone=2000
      maxiso=3000
   END
ENDCASE
RESOLVE_ROUTINE, $
   series+"nuc", $
   /COMPILE_FULL_FILE, $
   /NO_RECOMPILE  
CALL_PROCEDURE, $
   'get'+series+'runs', $
   mass, $
   METALLICITY=metallicity, $
   EXPLOSIONS=exp0

IF N_ELEMENTS(exp) EQ 0 THEN exp=exp0

dump="#nucleo"

nmass=N_ELEMENTS(mass)
nexp=N_ELEMENTS(exp)


FOR ie=0,nexp-1 DO BEGIN
; assure that these are zeroed
    xion=DBLARR(maxiso,maxzone+1,nmass)
    xmass=DBLARR(maxzone+2,nmass)
    xnzone=LONARR(nmass)
    xidx=LONARR(maxiso)

    FOR im=0,nmass-1 DO BEGIN
        ms=mass_string(mass[im])
        name=series+ms
        dir=dir0+name+'/'
        file=dir+name+exp[ie]+dump
        loaddump,file

        getdump,'nq 53',imaxb
        getdump,'nq 2',jm
        getdump,'xm',xm,CENTER=-1

        xnzone[im]=jm
      
; get zones - total mass in each ion  
        yion=DBLARR(imaxb,jm+1)
        FOR ii=0,jm-1 DO BEGIN
            getdump,'xb'+STRTRIM(STRING(ii+1),2),xb
            yion[*,ii]=xb*xm[ii]
        ENDFOR
        xmass[0:jm-1,im]=xm

 ; add wind
        getdump,'xwindb',xb
        yion[*,jm]=xb       
        getdump,'xmlost',xmwind
        xmass[jm,im]=xmwind
 
; piston mass
        getdump,'zm',zm,CENTER=-1
        xmass[jm+1,im]=zm[0]
       
; map ions
        getdump,'idxb',yidx
        niso=N_ELEMENTS(yidx)
        ymap=LONARR(niso)
        IF im EQ 0 THEN BEGIN
            ymap=INDGEN(niso)
            nxidx=niso
            xidx[0:niso-1]=yidx
        ENDIF ELSE BEGIN
            FOR ii=0,niso-1 DO BEGIN
                j=(WHERE(xidx[0:nxidx-1] EQ yidx[ii]))[0]
                IF j EQ -1 THEN BEGIN
                    xidx[nxidx]=yidx[ii]
                    ymap[ii]=nxidx
                    nxidx+=1
                ENDIF ELSE BEGIN
                    ymap[ii]=j
                ENDELSE
            ENDFOR            
        ENDELSE
        FOR ii=0,niso-1 DO BEGIN
            xion[ymap[ii],0:jm,im]=yion[ii,0:jm]
        ENDFOR
    ENDFOR

    print,ftag+'sorting-1'
    ii=SORT(xidx[0:nxidx-1])
    idxion=TEMPORARY(xidx[ii])
    xmaxzone=MAX(xnzone)
    print,ftag+'sorting-2'
    FOR iz=0,xmaxzone DO BEGIN
        xion[0:nxidx-1,iz,*]=xion[ii,iz,*]
    ENDFOR
    ;; xion =TEMPORARY( xion[ii,0:xmaxzone,*])

    print,ftag+'truncating'    
    xion =TEMPORARY( xion[0:nxidx-1,0:xmaxzone  ,*])
    xmass=TEMPORARY(xmass[          0:xmaxzone+1,*])

    nmass=LONG64(nmass)
    xmaxzone=LONG64(xmaxzone)
    maxalliso=LONG64(nxidx)
    idxion=LONG64(idxion)

    print,ftag+'writing' 
    version=10000LL
    OPENW, unit, dir0+series+'nuc'+exp[ie]+'.nucleo.zone.alliso.bgal',/GET_LUN,/SWAP_IF_BIG_ENDIAN
    WRITEU,unit,version
    WRITEU,unit,nmass,xmaxzone,maxalliso
    WRITEU,unit,idxion,xmass,xion
    CLOSE,unit
    FREE_LUN, unit
    print,ftag+'clearing' 
    clear,xion
    print,ftag+'Done.' 

ENDFOR

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO load_nucleo_zone_alliso_bgal,filename, $
                                 DIMENSIONS_ONLY=dimensions_only, $
                                 NMASS=nmass, $
                                 MAXZONE=maxzone, $
                                 MAXALLISO=maxalliso, $
                                 IDXION=idxion, $
                                 XMASS=xmass, $
                                 XION=xion                  

@compile_opt

starttime=SYSTIME(/SECONDS)

ftag=' [LOAD_NUCLEO_ZONE_ALLISO_BGAL] '

PRINT,ftag+'loading '+filename

OPENR, unit, filename, /GET_LUN, /SWAP_IF_BIG_ENDIAN
version=0LL
READU, unit,version
nmass=0LL
maxzone=0LL
maxalliso=0LL
READU,unit,nmass,maxzone,maxalliso

PRINT,ftag+'Version =   ',version
PRINT,ftag+'nmass =     ',nmass
PRINT,ftag+'maxzone =   ',maxzone
PRINT,ftag+'maxalliso = ',maxalliso

idxion=LON64ARR(maxalliso)
xmass=DBLARR(maxzone+2,nmass)
READU,unit,idxion,xmass

IF N_ELEMENTS(dimensions_only) EQ 0 THEN dimensions_only=0
IF dimensions_only EQ 0 THEN BEGIN 
    xion=DBLARR(maxalliso,maxzone+1,nmass)
    READU,unit,xion
ENDIF
CLOSE,unit
FREE_LUN, unit

finishtime=SYSTIME(/SECONDS)
PRINT,ftag+'data loaded in '+time2human(finishtime-starttime)

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO mix, yin, yout, xm, mbox0, iteration

@compile_opt

ftag=' [MIX] '

COMMON mix_save, $
  xm_save, $
  nzone_save, $
  mbox_save, $
  iteration_save, $
  mixmat_save

; get dimensions
yin=REFORM(temporary(yin))
dimensions=SIZE(yin,/DIMENSIONS)
nspec=dimensions[0]
nzone=dimensions[1]
IF N_ELEMENTS(dimensions) GT 2 THEN BEGIN
    nmodel=dimensions[2]
ENDIF ELSE BEGIN
    nmodel=1
ENDELSE

ii=WHERE(xm EQ 0)
IF ii[0] NE -1 THEN nzone=ii[0]

mbox=mbox0
@physconst
IF mbox LT 1.D20 THEN mbox*=xmsun

; remove remnant and wind
nzone-=2

new=0
IF N_ELEMENTS(xm_save) EQ 0 THEN BEGIN
    new=1
ENDIF ELSE IF N_ELEMENTS(xm_save) NE N_ELEMENTS(xm) THEN BEGIN 
    new=1
ENDIF ELSE IF nzone_save NE nzone THEN BEGIN
    new=1
ENDIF ELSE IF iteration_save NE iteration THEN BEGIN
    new=1
ENDIF ELSE IF TOTAL(((xm[0:nzone+1]-xm_save[0:nzone+1])/(xm[0:nzone+1]>1.D-99))^2) GT 1.D-10 THEN BEGIN
    new=1
ENDIF

IF new THEN begin
    mat=DBLARR(nzone,nzone)
; construct mixing matrix
    i=0
    xmt=0.0D0
    FOR j=i,nzone-1 DO BEGIN
        xmt+=xm[j]
        IF xmt GE mbox THEN BREAK
    ENDFOR
    IF xmt LT mbox THEN BEGIN
; complete mixing
        mat[0,0:nzone-1]=xm[0:nzone-1]/xmt
        GOTO,doneloop
    ENDIF
    mat[i,i:j]=xm[i]/xmt
    
    FOR i=1,nzone-1 DO BEGIN
        j0=j
        xmt-=xm[i-1]
        xmt0=xmt
        WHILE xmt LE mbox DO BEGIN
            j+=1
            xmt+=xm[j]
            IF j EQ nzone-1 THEN BREAK
        ENDWHILE
        IF i LE j0 THEN mat[i,0:j0]=mat[i-1,0:j0]*xmt0/xmt*xm[i]/xm[i-1]
        IF j GT j0 THEN mat[i,j0+1:j]=xm[i]/xmt
        
        IF j EQ nzone-1 THEN BREAK
    ENDFOR
doneloop:
    FOR ii=i+1,nzone-1 DO BEGIN
        mat[ii,0:nzone-1]=mat[i,0:nzone-1]*xm[ii]/xm[i]
    ENDFOR

; do mixing iterations
    IF N_ELEMENTS(iteration) EQ 0 THEN iteration=4
    mixmat=MATRIX_POWER(mat,iteration) 

; save new values
  xm_save=xm
  nzone_save=nzone
  mbox_save=mbox
  iteration_save=iteration
  mixmat_save=mixmat
ENDIF ELSE BEGIN
    PRINT,ftag+'Speedup: recycling mixing matrix.'
    mixmat=mixmat_save
ENDELSE

;  print,'he test'
;  x=yin[1,0:nzone-1]
;  y=TRANSPOSE(mixmat)##x
;  print,total(x)
;  print,total(y)

;  print,'fe test'
;  x=yin[25,0:nzone-1]
;  y=TRANSPOSE(mixmat)##x
;  print,total(x)
;  print,total(y)

;  print,'boundary test'
;  x=REPLICATE(0.D0,nzone)
;  x[nzone-1]=1.D0
;  y=TRANSPOSE(mixmat)##x
;  print,total(x)
;  print,total(y)

;  print,'boundary test'
;  FOR i=0,nzone-1 DO BEGIN
;      x=REPLICATE(0.D0,nzone)
;      x[i]=1.D0
;      y=TRANSPOSE(mixmat)##x
;      x=TOTAL(x)
;      y=TOTAL(y)
;      print,i,x,y,x-y
;  ENDFOR

;  tvscl,mixmat

yout=DBLARR(nspec,nzone+1,nmodel)
FOR i=0,nmodel-1 DO BEGIN
    IF !CPU.HW_NCPU EQ 1 THEN BEGIN
        yout[*,0:nzone-1,i]=TRANSPOSE(mixmat)##yin[*,0:nzone-1,i]
    ENDIF ELSE BEGIN
        yout[*,0:nzone-1,i]=MATRIX_MULTIPLY(yin[*,0:nzone-1,i],mixmat,/BTRANSPOSE)
    ENDELSE
    yout[*,nzone,i]=yin[*,nzone,i]
ENDFOR

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO masscut, yin, yout, xm, mcut0

@compile_opt

; get dimensions
yin=REFORM(temporary(yin))
dimensions=SIZE(yin,/DIMENSIONS)
nspec=dimensions[0]
nzone=dimensions[1]
IF N_ELEMENTS(dimensions) GT 2 THEN BEGIN
    nmodel=dimensions[2]
ENDIF ELSE BEGIN
    nmodel=1
ENDELSE

ii=WHERE(xm EQ 0)
IF ii[0] NE -1 THEN nzone=ii[0]

; remove remnant (but not wind)
nzone-=1

@physconst
mcut=mcut0
IF mcut LT 1.D20 THEN mcut*=xmsun

map=DBLARR(nzone)
xmt=xm[nzone]
i=0
WHILE xmt LT mcut DO BEGIN
    xmt+=xm[i]
    i+=1
ENDWHILE
i-=1
IF i GE 0 THEN map[i]=(xmt-mcut)/xm[i]
map[i+1:nzone-1]=1.D0

; print,i,map[i],mcut

yout=DBLARR(nspec,nmodel)
FOR i=0,nmodel-1 DO BEGIN
    yout[*,i]=map##yin[*,0:nzone-1,i]
ENDFOR

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO get_he_core_size, HECORE=hecore,SERIES=series,MASS=mass

@compile_opt

IF N_ELEMENTS(series) EQ 0 THEN series='u'
file='/u/alex/kepler/'+series+'nuc/'+series+'nuc.gal'
loadgal,file
common gal, ngalidx,ngalmass,galidx,galmass,galyield,remgal
common gal3, $
  begal, begalname, $
  mcoregal, mcoregalname, $
  rcoregal, rcoregalname, $1
  tcoregal, tcoregalname, $
  dcoregal, dcoregalname, $
  scoregal, scoregalname

ii=WHERE(mcoregalname EQ 'He core')
ii=ii[0]
hecore=mcoregal[*,ii]
mass=galmass
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION energy_interpolate, data, energy_set, mass_set, parameters

@compile_opt

; defaults
energy=1.2D51
mode=0
exponent=0.0D0

; parameters can be array with values in order of above list
; or a structure with tags 
;{ENERGY:1.2D51, MODE:0, EXPONENT:0.D0}
IF N_ELEMENTS(parameters) NE 0 THEN BEGIN
    IF SIZE(parameters,/TYPE) NE 8 THEN BEGIN
        n=N_ELEMENTS(parameters)<3
        SWITCH n OF
            3:exponent=parameters[2]
            2:mode=LONG(parameters[1])
            1:energy=parameters[0]
        ENDSWITCH
    ENDIF ELSE BEGIN
        n=N_TAGS(parameters)
        names=TAG_NAMES(parameters)
        FOR i=0,n-1 DO BEGIN
            x=parameters.(i)
            CASE names[i] OF
                'ENERGY': energy=x
                'MODE': mode=LONG(x)
                'EXPONENT': exponent=x
            ENDCASE
        ENDFOR
    ENDELSE
ENDIF

IF N_ELEMENTS(mode) EQ 0 THEN mode=0
CASE mode of 
    0:BEGIN
        chi=((energy_set-energy)/energy)^2
        chi_x=SORT(chi)
        i0=chi_x[0]
        i1=chi_x[1]
        e0=energy_set[i0]
        e1=energy_set[i1]
        de=e1-e0
        d0=energy-e0
        d1=e1-energy
        f0=d1/de
        f1=d0/de
        IF chi[0] LT 1.d-3 THEN BEGIN
            f0=1.D0
            f1=0.D0
        ENDIF
        IF (f0 LT 0) OR (f1 GT 1) THEN BEGIN
            f0=0
            f1=1
        ENDIF
        IF (f1 LT 0) OR (f0 GT 1) THEN BEGIN
            f0=1
            f1=0
        ENDIF
        x=size(data,/DIMENSIONS)
        n=PRODUCT(x[0:N_ELEMENTS(x)-2],/INTEGER)
        y=[n,x[N_ELEMENTS(x)-1]]        
        val=REFORM(data,y)
        val=val[*,i0]*f0+val[*,i1]*f1
        y=x[0:N_ELEMENTS(x)-2]
        val=REFORM(val,y)
    END
    1:BEGIN
        e=energy*(mass_set/20.D0)^exponent
        n=N_ELEMENTS(mass_set)
        val=data[*,*,0]
        FOR i=0,n-1 DO BEGIN
            chi=((energy_set-e[i])/e[i])^2
            chi_x=SORT(chi)
            i0=chi_x[0]
            i1=chi_x[1]
            e0=energy_set[i0]
            e1=energy_set[i1]
            de=e1-e0
            d0=e[i]-e0
            d1=e1-e[i]
            f0=d1/de
            f1=d0/de
            IF (f0 LT 0) OR (f1 GT 1) THEN BEGIN
                f0=0
                f1=1
            ENDIF
            IF (f1 LT 0) OR (f0 GT 1) THEN BEGIN
                f0=1
                f1=0
            ENDIF
            val[*,i]=data[*,i,i0]*f0+data[*,i,i1]*f1
        ENDFOR
    END
ENDCASE

RETURN, val

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION yield_interpolate, data, coordinate, value

@compile_opt

ftag=' [YIELD_INTERPOLATE] '

x=size(data,/DIMENSIONS)
n=PRODUCT(x[0:N_ELEMENTS(x)-2],/integer)
y=[n,x[N_ELEMENTS(x)-1]]

val=REFORM(data,y)

chi=(coordinate-value)^2
chi_x=SORT(chi)
i0=chi_x[0]
i1=chi_x[1]
e0=coordinate[i0]
e1=coordinate[i1]
de=e1-e0
d0=value-e0
d1=e1-value
f0=d1/de
f1=d0/de
IF chi[0] LT 1.d-3 THEN BEGIN
    f0=1.D0
    f1=0.D0
ENDIF

eps=1.D-7
IF (f1 GT 1.D0+eps) or $
  (f1 LT 0.D0-eps) OR $
  (f0 GT 1.D0+eps) OR $
  (f0 LT 0.D0-eps) THEN BEGIN
    PRINT,ftag+'WARNING: OUT OF RANGE! - TRUNCATING'
ENDIF

f0>=0.D0
f0<=1.D0
f1>=0.D0
f1<=1.D0
val=val[*,i0]*f0+val[*,i1]*f1

y=x[0:N_ELEMENTS(x)-2]
val=REFORM(val,y)

RETURN,val
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO obsolete_load_mixlib,filename, $
                NEXP=nexp, $
                NMASS=nmass, $
                NMIX=nmix, $
                NISO=niso, $
                IDXION=idxion, $
                ENERGY=energy, $
                MIXBOX=mixbox, $
                MASS=mass, $
                XISOM=xisom

@compile_opt

ftag=' [LOAD_MIXLIB] '

starttime=SYSTIME(/SECONDS)
PRINT,ftag+'loading '+filename

OPENR, unit, filename, /GET_LUN,/SWAP_IF_BIG_ENDIAN
version=0LL
READU, unit,version
nmix=0LL
nexp=0LL
nmass=0LL
niso=0LL
READU,unit,niso,nmass,nexp,nmix

PRINT,ftag+'Version =   ',version
PRINT,ftag+'nexp =      ',nexp
PRINT,ftag+'nmass =     ',nmass
PRINT,ftag+'nmix =      ',nmix
PRINT,ftag+'niso =      ',niso

idxion=LON64ARR(niso)
mass=DBLARR(nmass)
energy=DBLARR(nexp)
mixbox=DBLARR(nmix)

xisom=DBLARR(niso,nmass,nexp,nmix)

READU,unit,idxion,mass,energy,mixbox
READU,unit,xisom
CLOSE,unit
FREE_LUN, unit

finishtime=SYSTIME(/SECONDS)
PRINT,ftag+'data loaded in '+time2human(finishtime-starttime)

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO yieldlib, $
   SERIES=series

IF N_ELEMENTS(series) EQ 0 THEN series='z'

;; getsize, $
;;    series
make_nucleo_zone_alliso_bgal, $
   SERIES=series

galdata,SERIES=series,/presn,/entropy,/radius,/density,/temperature,/core,/be
; DOES IT LOAD OTHER DATA AS WELL ????
; galdata,series='z',['A','B','C','D','E','F','G','H','I','J','P','V']

cuts=['Ye','S4']

;; FOR i=0,N_ELEMENTS(cuts)-1 DO BEGIN
;;    make_alliso_mixlib, SERIES=series, CUT=cuts[i]
;;    make_radiso_mixlib, SERIES=series, CUT=cuts[i]
;;    make_deciso_mixlib, SERIES=series, CUT=cuts[i]
;;    make_el_mixlib,     SERIES=series, CUT=cuts[i]
;; ENDFOR

makestardb,SERIES=series

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------