PRO maxradio, $
   ions, $                ; INPUT:  ion list
   map,$                  ; OUTPUT: map matrix
   radions, $             ; OUTPUT: ion list
   VERBOSE=verbose, $     ; write out more stuff
   SILENT=silent, $       ; do not write out anything
   STABLE=stable, $       ; only output stable ions, discrad decay chains
   IONLIST=ionlist, $     ; only output these isotops
   ELEMENTS=elements, $   ; map to elements
   MASS=mass, $           ; map to isobars
   SOLPROD=solprod, $     ; return solar prduction factor matrix
   KEEPALL=keepall, $     ; do not discard ions not needed, to ensure always the same complete ion list (e.g., solar stable isotopes) is returend 
   ABUNDANCE=abundance, $ ; map abundance (mol fractions) not mass fractions
   RADIONA=radiona, $     ; OUTPUT: return ion mass number A
   RADIONZ=radionz, $     ; OUTPUT: return ion charge number Z
   DECAYFILE=decayfile, $ ; overwrite default foe decay file
   X2Y=x2y, $             ; map from mass fraction input to mol fraction output
   ADDIONS=addions, $     ; INPUT: also output these extra ions, but not add to EL/A sums
   KEEPIONS=keepions, $   ; INPUT: ions to keep as stable but not add ion EL/A sums
   STABIONS=stabions      ; INPUT: ions to fully treat as stable
  

; the following is the final plan; this is not implemented as of 20100425.
;-----------------------------------------------------------------------
;            Treat as     Keep in    Add to      compatible
;            stable       output     EL/A sum    with solabu                 
;            (no decays)
; 
; ionlist    DEPENDS      YES        DEPENDS     CHECK 
; addions    NO           YES        NO          NO
; keepions   YES          YES        NO          NO
; stabions   YES          YES        YES         IONS: NO; EL/A:CHECK
;-----------------------------------------------------------------------

@compile_opt
ftag=' [MAXRADIO] '

IF N_ELEMENTS(x2y) EQ 0 THEN x2y=0
IF N_ELEMENTS(abundance) EQ 0 THEN abundance=0
IF N_ELEMENTS(verbose) EQ 0 THEN verbose=1
IF N_ELEMENTS(silent) EQ 0 THEN silent=0
IF silent EQ 1 THEN verbose=0

IF N_ELEMENTS(decayfile) EQ 0 THEN BEGIN
   decayfile='/u/alex/kepler/local_data/decay.dat'
ENDIF

decays=['s','a','b-','b+','n','p','a2','p2']
da    =[  0, -4,   0,   0, -1, -1,  -8,  -2]
dz    =[  0, -2,  +1,  -1,  0, -1,  -4,  -2]
arem  =[  0,  4,   0,   0,  1,  1,   4,   1]
zrem  =[  0,  2,   0,   0,  0,  1,   2,   1]
frem  =[  0,  1,   0,   0,  1,  1,   2,   2]

nd=N_ELEMENTS(decays)

ni=N_ELEMENTS(ions)
aions=LONARR(ni)
zions=LONARR(ni)
iidex=LONARR(ni)
unsorted=0

FOR i=0,ni-1 DO BEGIN
    ionaz,ions[i],aionx,zionx
    aions[i]=aionx
    zions[i]=zionx    
    iidex[i]=aionx+1000L*zionx
    IF i GT 0 THEN BEGIN
        IF (iidex[i-1] GT iidex[i]) THEN BEGIN
            IF verbose EQ 1 THEN PRINT,ftag+'WARNING: ion array not proporly sorted.'
            unsorted=1
        ENDIF
    ENDIF
ENDFOR
aionmax=MAX(aions)
zionmax=MAX(zions)

IF unsorted THEN BEGIN
    IF verbose EQ 1 THEN PRINT,ftag+'sorting ions...'
    sortmap=SORT(iidex)
    iidex=iidex[sortmap]
    aions=aions[sortmap]
    zions=zions[sortmap]
    sortmap=SORT(sortmap) ; make inverse sort map
ENDIF

mr=MAX([2000L,3L*ni])
ar=INTARR(mr)
zr=INTARR(mr)
br=DBLARR(mr,nd)
irdex=LONARR(mr)

;-----------------------------------------------------------------------
; read DECAY file
;-----------------------------------------------------------------------
OPENR,unit,decayfile,/GET_LUN

nr=0
WHILE NOT EOF(unit) DO BEGIN
    l=''
    READF,unit,l
    IF STRMID(l,0,1) EQ ';' THEN GOTO,end_of_loop
    a=STRSPLIT(l,/EXTRACT)
    el=a[0]
    ionaz,el,ax,zx
;    IF (ax GT aionmax) AND (zx GT zionmax) THEN GOTO,end_of_loop
    IF (ax GT aionmax) THEN GOTO,end_of_loop
    na=N_ELEMENTS(a)
    IF na EQ 2 THEN f=1.D0 ELSE f=DOUBLE(a[1])
    d=a[na-1]
    ix=(WHERE(d EQ decays))[0]
    IF ix EQ -1 THEN BEGIN
        PRINT,ftag+'unknown decay in '+l
    ENDIF
    xrdex=ax+1000L*zx
    ii=1
    IF nr GT 0 THEN BEGIN
        IF xrdex EQ irdex[nr-1] THEN ii=0 
    ENDIF
    IF ii EQ 1 THEN BEGIN
        ar[nr]=ax
        zr[nr]=zx
        irdex[nr]=ax+1000L*zx
        nr+=1
    ENDIF
    br[nr-1,ix]=f
end_of_loop:
ENDWHILE

CLOSE,unit
FREE_LUN,unit
;-----------------------------------------------------------------------

IF verbose NE 0 THEN print,ftag+'nr=',nr

; add all isotopes from burn
j=0
FOR i=0,ni-1 DO BEGIN
beg_loop_addburn:
    IF j LT nr THEN BEGIN
        IF irdex[j] LT iidex[i] THEN BEGIN
            j=j+1
            GOTO,beg_loop_addburn
        ENDIF
        IF irdex[j] EQ iidex[i] THEN GOTO,end_loop_addburn
    ENDIF
    IF j EQ nr THEN BEGIN
        jj=j-1
    ENDIF ELSE BEGIN
        ar[j+1:nr]=ar[j:nr-1]
        zr[j+1:nr]=zr[j:nr-1]
        br[j+1:nr,*]=br[j:nr-1,*]
        irdex[j+1:nr]=irdex[j:nr-1]        
        IF zr[j] EQ zions[i] THEN BEGIN
            jj=j
        ENDIF ELSE BEGIN
            jj=j-1
        ENDELSE
    ENDELSE
    ar[j]=aions[i]
    zr[j]=zions[i]
    irdex[j]=iidex[i]
    br[j,*]=br[jj,*]
    nr=nr+1
end_loop_addburn:
ENDFOR

IF verbose NE 0 THEN print,ftag+'nr=',nr

; add all decayed isotopes.
jj=0
REPEAT BEGIN
    inew=0
    FOR jj=0,nr-1 DO BEGIN
        FOR i=1,nd-1 DO BEGIN
            IF br[jj,i] NE 0.0D0 THEN BEGIN
                an=ar[jj]+da[i]
                zn=zr[jj]+dz[i]
                IF (an LE aionmax) AND (zn LE zionmax) AND (an GT 0) AND (zn GE 0) THEN BEGIN
                    in=an+1000L*zn
                    IF in GT irdex[j] THEN BEGIN
                        WHILE irdex[j] LT in DO j=j+1
                        IF irdex[j] EQ in THEN GOTO,end_loop_br
                        IF zr[j] EQ zr[jj] THEN BEGIN
                            ji=j+1
                        ENDIF ELSE BEGIN
                            ji=j-1
                        ENDELSE
                        k=0
                    ENDIF ELSE BEGIN
                        WHILE irdex[j] GT in DO j=j-1
                        IF irdex[j] EQ in THEN GOTO,end_loop_br
                        IF zr[j] EQ zr[jj] THEN BEGIN
                            ji=j
                        ENDIF ELSE BEGIN
                            ji=j+2
                        ENDELSE
                        k=1
                    ENDELSE
                    ar[j+1+k:nr]=ar[j+k:nr-1]
                    zr[j+1+k:nr]=zr[j+k:nr-1]
                    br[j+1+k:nr,*]=br[j+k:nr-1,*]
                    irdex[j+1+k:nr]=irdex[j+k:nr-1]        
                    ar[j+k]=an
                    zr[j+k]=zn
                    irdex[j+k]=in
                    br[j+k,*]=br[ji,*]
                    nr=nr+1
                    inew=inew+1                
                ENDIF
            ENDIF
end_loop_br:
        ENDFOR
    ENDFOR
IF verbose NE 0 THEN PRINT,ftag+'loop ',inew
ENDREP UNTIL inew EQ 0

IF verbose NE 0 THEN PRINT,ftag+'nr=',nr

;FOR i=0,nr-1 DO BEGIN
;    s=''
;    FOR j=0,nd-1 DO BEGIN
;        IF br(i,j) NE 0.0D0 THEN s=s+decays(j)
;    ENDFOR
;    PRINT,irdex(i),ar(i),zr(i),s
;ENDFOR

; indices of decay products
tr=REPLICATE(-1L,nr,nd)
FOR i=0,nr-1 DO BEGIN
    FOR j=0,nd-1 DO BEGIN
        IF (j EQ 0) OR (br[i,j] NE 0.0D0) THEN BEGIN
            an=ar[i]+da[j]
            zn=zr[i]+dz[j]
            IF (j NE 0) AND (an EQ 0) AND (zn EQ 0) THEN PRINT,ftag+'error0'
            IF (an LE aionmax) AND (zn LE zionmax) THEN BEGIN
                in=an+1000L*zn
                jj=i
                IF in GT irdex[jj] THEN BEGIN
                    WHILE irdex[jj] LT in DO jj=jj+1
                    IF irdex[jj] NE in THEN PRINT,ftag+'error1'
                ENDIF ELSE BEGIN
                    WHILE irdex[jj] GT in DO jj=jj-1
                    IF irdex[jj] NE in THEN PRINT,ftag+'error2'
                ENDELSE
                tr[i,j]=jj
            ENDIF
        ENDIF
    ENDFOR
ENDFOR

IF verbose NE 0 THEN PRINT,ftag+'done index1'

sr=LONARR(ni)
; indices of original ions into decay ions
j=0
FOR i=0,ni-1 DO BEGIN
beg_loop_ind:
    IF irdex[j] LT iidex[i] THEN BEGIN
        j=j+1
        GOTO,beg_loop_ind
    ENDIF
    IF irdex[j] NE iidex[i] THEN  BEGIN
        PRINT,ftag+'error3 (probably He4 is not properly sorted in)'
    END
    sr[i]=j
ENDFOR

IF verbose NE 0 THEN PRINT,ftag+'done index2'

irem=REPLICATE(-1,nd)
FOR i=1,nd-1 DO BEGIN
    IF da[i] NE 0 THEN BEGIN
        ii=arem[i]+1000L*zrem[i]
        FOR j=0,nr-1 DO BEGIN
            IF irdex[j] EQ ii THEN BEGIN
                irem[i]=j
                GOTO,end_loop_irem
            ENDIF 
        ENDFOR
    ENDIF
end_loop_irem:
ENDFOR

IF verbose NE 0 THEN PRINT,ftag+'done index3'

; make dec matrix
dec=DBLARR(ni,nr)
stack=INTARR(nr)
fstack=DBLARR(nr)

FOR i=0,ni-1 DO BEGIN
    istack=1
    fstack[0]=1.0D0
    stack[0]=sr[i]

    WHILE (istack NE 0) DO BEGIN
        istack=istack-1
        item=stack[istack]
        ftem=fstack[istack]
        IF abundance EQ 0 THEN ftem=ftem*ar[item]/(1.D0*aions[i])
        dec[i,item]=dec[i,item]+ftem

        FOR j=1,nd-1 DO BEGIN
            IF tr[item,j] GE 0 THEN BEGIN
                trx=tr[item,j]
                brx=br[item,j]
                stack[istack]=trx
                fstack[istack]=brx
                istack=istack+1
                IF irem[j] NE -1 THEN BEGIN
                    trx=irem[j]
                    stack[istack]=trx
                    fstack[istack]=brx*frem[j]                    
                    istack=istack+1
                ENDIF
            ENDIF
        ENDFOR        
    ENDWHILE
ENDFOR

IF verbose NE 0 THEN PRINT,ftag+'done dec matrix'

IF N_ELEMENTS(stable) EQ 0 THEN BEGIN 
    IF N_ELEMENTS(addions) NE 0 THEN stable=1 ELSE stable=0
ENDIF
IF N_ELEMENTS(solprod) EQ 0 THEN solprod=0
IF N_ELEMENTS(elements) EQ 0 THEN elements=0
IF elements NE 0 THEN BEGIN
    stable=1   
ENDIF
IF N_ELEMENTS(mass) EQ 0 THEN mass=0
IF mass NE 0 THEN BEGIN
    stable=1   
ENDIF
IF (mass NE 0) AND (elements NE 0) THEN BEGIN
    PRINT,ftag+'ERROR: "mass" and "elements" specified...skipping mass.'
    mass=0
ENDIF
IF (stable NE 0) OR (N_ELEMENTS(ionlist) NE 0) THEN BEGIN
    @solabu_var
    IF (N_ELEMENTS(ionlist) EQ 0) OR (elements NE 0) THEN BEGIN
        IF N_ELEMENTS(nsoldata) EQ 0 THEN solabu
        ionlist=sol_ion
        IF (N_ELEMENTS(addions) NE 0) THEN BEGIN
            IF (elements EQ 0) AND (mass EQ 0) THEN BEGIN
                ionlist=[ionlist,addions]
            ENDIF ELSE BEGIN
                PRINT,ftag+'ERROR: "addions" specified...skipping "elements/mass"'
            ENDELSE
        ENDIF
    ENDIF
    ionaz,ionlist,outa,outz
    iiout=1000L*outz+outa
    iiout=iiout[UNIQ(iiout,SORT(iiout))]
    ns=N_ELEMENTS(iiout)
    outmap=REPLICATE(-1,ns)    
    j=0
    i=0
    WHILE (j LT ns) AND (i LT nr) DO BEGIN
       WHILE (iiout[j] LT irdex[i]) DO BEGIN
          outmap[j]=-1
          j+=1
          if j GE ns THEN GOTO, end_map
       ENDWHILE     
       IF (iiout[j] EQ irdex[i]) THEN BEGIN
          outmap[j]=i
          j+=1
       ENDIF
       i+=1
    ENDWHILE
    end_map:

    IF N_ELEMENTS(keepall) EQ 0 THEN keepall=0
    IF keepall EQ 1 THEN BEGIN
        x=WHERE(outmap EQ -1)
        IF x[0] NE -1 THEN BEGIN
            outmap[x]=nr
            dec=[[TEMPORARY(dec)],[REPLICATE(0.0D0,ni)]]
        ENDIF
    ENDIF

    x=WHERE(outmap NE -1)
; here we still have to deal with the case that there is no output,
; all ions removed ... Often we may want to use /KEEPALL
    IF x[0] EQ -1 THEN BEGIN
       PRINT,ftag+'ERROR: no isotopes to output.'
       RETURN
    ENDIF
    outmap=outmap[x]
    dec=dec[*,outmap]
    IF verbose NE 0 THEN PRINT,ftag+'done stable dec matrix'
    ar=iiout mod 1000L
    zr=iiout  /  1000L
    nr=N_ELEMENTS(x)


    IF (elements NE 0) THEN BEGIN
        m=REPLICATE(0L,1000)
        j=-1
        en0=-1
        FOR i=0,nr-1 DO BEGIN
            en1=zr[i]
            IF en1 NE en0 THEN BEGIN
                j+=1    
                dec[*,j]=dec[*,i]
            ENDIF ELSE BEGIN
                dec[*,j]+=dec[*,i]
            ENDELSE
            en0=en1
            m[j]=i
        ENDFOR
        m=m[0:j]
        nr=j+1
        zr=zr[m]
        ar=zr*0-1
        dec=TEMPORARY(dec[*,0:j])

        IF (solprod NE 0) THEN BEGIN
            IF N_ELEMENTS(nsoldata) EQ 0 THEN solabu
            solidx=sol_el_z
            j=0
            FOR i=0,nr-1 DO BEGIN
                WHILE (solidx[j] LT zr[i]) AND (j LT (nsol_el-2)) DO j=j+1
                IF solidx[j] EQ zr[i] THEN BEGIN
                    IF abundance NE 0 THEN BEGIN 
                        fsol=1.D0/sol_el_mol[j]
                    ENDIF ELSE BEGIN
                        fsol=1.D0/sol_el_abu[j]
                    ENDELSE
                    dec[*,i]=dec[*,i]*fsol
                ENDIF ELSE BEGIN
                    PRINT,ftag+'ERROR: no solar abundance for element ',zr[i]
                    dec[*,i]=dec[*,i]*0.0D0
                ENDELSE
            ENDFOR
         ENDIF
    ENDIF ELSE IF (mass NE 0) THEN BEGIN
        m=REPLICATE(-1L,nr)
        FOR i=0,nr-1 DO BEGIN
            j=ar[i]
            IF (m[j] EQ -1) THEN BEGIN
                m[j]=i
            ENDIF ELSE BEGIN
                jj=m[j]
                dec[*,jj]=dec[*,jj]+dec[*,i]
            ENDELSE
        ENDFOR
        x=WHERE(m NE -1)
        m=m[x]
        ar=ar[m]
        dec=TEMPORARY(dec[*,m])
        zr=REPLICATE(0,nr)    
        nr=N_ELEMENTS(m)

        IF (solprod NE 0) THEN BEGIN
           IF N_ELEMENTS(nsoldata) EQ 0 THEN solabu
           solidx=sol_aa_a
           j=0 
           FOR i=0,nr-1 DO BEGIN
              WHILE (solidx[j] LT ar[i]) AND (j LT (nsol_aa-2)) DO j=j+1
              IF solidx[j] EQ ar[i] THEN BEGIN
                 IF abundance NE 0 THEN BEGIN 
                    fsol=1.D0/sol_aa_mol[j]
                 ENDIF ELSE BEGIN
                    fsol=1.D0/sol_aa_abu[j]
                 ENDELSE
                 dec[*,i]=dec[*,i]*fsol
              ENDIF ELSE BEGIN
                 PRINT,ftag+'ERROR: no solar abundance for mass ',ar[i]
                 dec[*,i]=dec[*,i]*0.0D0
              ENDELSE
           ENDFOR
        ENDIF 
     ENDIF ELSE BEGIN
        IF (solprod NE 0) THEN BEGIN
            IF N_ELEMENTS(nsoldata) EQ 0 THEN solabu
            solidx=sol_a+1000*sol_z
            j=0
            FOR i=0,nr-1 DO BEGIN
                WHILE (solidx[j] LT iiout[i]) AND (j LT (nsoldata-2)) DO j=j+1
                IF solidx[j] EQ iiout[i] THEN BEGIN
                    IF abundance NE 0 THEN BEGIN 
                        fsol=1.D0/sol_mol[j]
                    ENDIF ELSE BEGIN
                        fsol=1.D0/sol_abu[j]
                    ENDELSE
                    dec[*,i]=dec[*,i]*fsol
                ENDIF ELSE BEGIN
                    PRINT,ftag+'ERROR: no solar abundance for ',iiout[i]
                    dec[*,i]=dec[*,i]*0.0D0
                ENDELSE
            ENDFOR
         ENDIF 
    ENDELSE

ENDIF

; convert mass fraction of isotopes to element mole fractions
; NOTE: for just mapping of isotope mole fractions to element mole
; fractions do not use the ABUNDANCE keyword; the result from mass
; fraction is fine.
IF (solprod EQ 0) AND x2y NE 0 THEN diag_mult,dec,1.D0/aions,1

map=TEMPORARY(dec)
radiona=ar[0:nr-1]
radionz=zr[0:nr-1]
radions=ionstr(radiona,radionz)

IF unsorted THEN BEGIN
    map=map[sortmap,*]
ENDIF

END