; helper function
FUNCTION spaghetti_nolabel, a1,a2,a3,a4
RETURN,""
END


PRO spaghetti,x,y,ny, $
              XMULTI=xmulti, $ ; x is also a multi-D array
              NN=nn, $ ; array of max index values
              XXFRAC=xxfrac, $ ; fractional (incremental) locations of breaks on the x-axis
              XXRANGE=xxrange, $ ; incremetal xrange stops
              XXTICKLEN=xxticklen, $ ; length of the tick mark drawn on the break -- major tickmark if omitted or zero; 1 draws full y-axis; default !P.TICKLEN
              XXANO=xxano, $ ; annotation writen below the break
              XXNOTICK=xxnotick, $ ; suppess tickmarks on break
              YLOG=ylog, $ 
              YRANGE=yrange, $
              XTITLE=xtitle, $
              YTITLE=ytitle, $
              XMARGIN=xmargin, $
              YMARGIN=ymargin, $
              LINE=line, $ ; array of line styles
              COLOR=color, $ ; array of colors
              THICK=thick, $ ; array of line thicknesses
              XXTICKS=xxticks, $ ; number of tickmark subsections (makes n+1 trickmarks) [usually not useful]
              XXTICKINTERVAL=xxtickinterval, $ ; size of interval used for tickmarks
              XXLABEL=xxlabel, $ ; label to be drawn instead of tickmarks
              XXTITLE=xxtitle0, $ ; title to be drawn under each segment
              XXTICKFORMAT=xxtickformat, $ ; xtickformat for indiviual sections
              ONELINE=oneline, $ ; draw contineous lines rather than separate lines for each section
              PROC=proc, $
              XMAP=xmap, $ ; need to set /oneline to get mapping data
              YMAP=ymap, $ ; use /nodata to suppress plotting
              NOERASE=noerase, $
              NODATA=nodata, $
              NOFRAME=noframe, $
              REDUCE=reduce, $
              SETNORMAL=setnormal, $ ; do switch to normal coordinate system 
              NONORMAL=nonarmal, $ ; do not switch to normal coordinate system 
              NORMALONLY=normalonly, $ ; do only switch to normal coordinate system 
              BACKGROUNDCOLOR=backgroundcolor, $
              OUTSIDECOLOR=outsidecolor, $
              TITLE=title, $
              XTANGLE=xtangle, $
              LEFTAXIS=leftaxis, $
              RIGHTAXIS=rightaxis, $
              XXEDGECOLOR=xxedgecolor



; x=indgen(1001)*0.02
; y=0.5*(1+0.997*Sin((x*10)^0.5))

IF N_ELEMENTS(ny) EQ 0 THEN BEGIN
    nydn=SIZE(y,/N_DIMENSIONS)
    nyd=SIZE(y,/DIMENSIONS)
    IF nydn EQ 1 THEN ny=1 ELSE ny=nyd[1]
ENDIF 
IF N_ELEMENTS(xmulti) EQ 0 THEN BEGIN
    nxdn=SIZE(x,/N_DIMENSIONS)
    nxd=SIZE(x,/DIMENSIONS)
    IF nxdn EQ 1 THEN BEGIN
        nx=1 
        xmulti=0
    ENDIF ELSE BEGIN
        nx=nxd[1]
        IF nx NE ny THEN BEGIN
            PRINT,' [SPAGHETTY] ERROR: nx='+STRING(nx)+', ny='+STRING(ny)
            RETURN
        ENDIF
        xmulti=1
    ENDELSE
ENDIF

IF N_ELEMENTS(NN) EQ 0 THEN BEGIN
    nydn=SIZE(y,/N_DIMENSIONS)
    nyd=SIZE(y,/DIMENSIONS)
    IF nydn EQ 1 THEN nyn=1 ELSE nyn=nyd[1]
    nn=REPLICATE(nyd[0],nyn)
ENDIF

IF N_ELEMENTS(xtitle) EQ 0 THEN xtitle='interior mass / solar masses'
IF N_ELEMENTS(ytitle) EQ 0 THEN ytitle='mass fraction'

IF N_ELEMENTS(ylog) EQ 0 THEN ylog=0

IF N_ELEMENTS(xxfrac) EQ 0 THEN xxfrac=[0,0.4,0.8,1.0]

n=N_ELEMENTS(xxfrac)-1

IF N_ELEMENTS(xxrange) EQ 0 THEN xxrange=[0,1.6,6,20]
IF N_ELEMENTS(xxticks) EQ 0 THEN xxticks=0
IF N_ELEMENTS(xxticks) LT n THEN xxticks=REPLICATE(xxticks[0],n)

IF N_ELEMENTS(yrange) EQ 0 THEN BEGIN 
    ymin=MIN(y,MAX=ymax)
    yrange=[ymin,ymax]
ENDIF

IF N_ELEMENTS(noerase) EQ 0 THEN noerase=0

IF N_ELEMENTS(oneline) EQ 0 THEN oneline=0
IF N_ELEMENTS(nodata) EQ 0 THEN nodata=0
IF N_ELEMENTS(noframe) EQ 0 THEN noframe=0

IF N_ELEMENTS(reduce) EQ 0 THEN reduce=1.D-4

IF oneline EQ 1 THEN BEGIN
    y1=y*1.D0
    x1=x*1.D0
    last1=REPLICATE(-1L,nx)
    end1=REPLICATE(-1L,nx)
ENDIF

IF N_ELEMENTS(color) EQ 0 THEN color=!P.COLOR
IF N_ELEMENTS(color) LT ny THEN color=REPLICATE(color[0],ny)
IF N_ELEMENTS(line) EQ 0 THEN line=0
IF N_ELEMENTS(line) LT ny THEN line=REPLICATE(line[0],ny)
IF N_ELEMENTS(thick) EQ 0 THEN thick=!P.THICK
IF N_ELEMENTS(thick) LT ny THEN thick=REPLICATE(thick[0],ny)

IF N_ELEMENTS(xxticklen) EQ 0 THEN xxticklen=!P.TICKLEN
IF N_ELEMENTS(xxticklen) LT n-1 THEN xxticklen=REPLICATE(xxticklen,n-1) 

IF N_ELEMENTS(xxnotick) EQ 0 THEN xxnotick=0
IF N_ELEMENTS(xxnotick) LT n-1 THEN xxnotick=REPLICATE(xxnotick,n-1) 

;IF N_ELEMENTS() EQ 0   THEN =0
;IF N_ELEMENTS() LT n-1 THEN =REPLICATE(,n-1) 
IF N_ELEMENTS(xxtickinterval) EQ 0   THEN xxtickinterval=0
IF N_ELEMENTS(xxtickinterval) LT n THEN xxtickinterval=REPLICATE(xxtickinterval,n) 

IF N_ELEMENTS(xxtickformat) EQ 0 THEN xxtickformat=''
IF N_ELEMENTS(xxtickformat) LT n THEN xxtickformat=REPLICATE(xxtickformat,n) 

ymargin0=0
IF xtitle NE '' AND N_ELEMENTS(xxtitle0) NE 0 THEN ymargin0=1.5

IF N_ELEMENTS(xxtitle0) EQ 0 THEN xxtitle='' ELSE xxtitle=xxtitle0
IF N_ELEMENTS(xxtitle) LT n THEN xxtitle=REPLICATE(xxtitle,n) 

IF N_ELEMENTS(xmargin) EQ 0 THEN xmargin=[8.3,0]
IF N_ELEMENTS(ymargin) EQ 0 THEN ymargin=[3.1+ymargin0,0.5*(SQRT(5.)-1)]

chsize=CONVERT_COORD(!D.X_CH_SIZE,!D.Y_CH_SIZE,/DEVICE,/TO_NORMAL)
chsize=chsize(0:1)

frame0=[0,0,1,1]
frame=frame0+[+xmargin(0),+ymargin(0),-xmargin(1),-ymargin(1)]*chsize([0,1,0,1])

width=frame(2)-frame(0)

delta=chsize(0)
bheight=chsize(1)*0.5*(SQRT(5.)-1)

IF N_ELEMENTS(backgroundcolor) EQ 0 THEN backgroundcolor=!P.BACKGROUND
IF N_ELEMENTS(outsidecolor) EQ 0 THEN outsidecolor=!P.BACKGROUND
IF N_ELEMENTS(xxedgecolor) LT n-1 THEN xxedgecolor=REPLICATE(backgroundcolor,n-1)


IF noerase EQ 0 THEN BEGIN
    plot,[0,1],/NODATA, $
      POSITION=frame,/NORMAL, $
      YSTYLE=5, $
      XSTYLE=5 
    IF N_ELEMENTS(fill) EQ 0 THEN fill=1
    IF fill EQ 1 THEN BEGIN
        POLYFILL,[0,0,1,1],[0,1,1,0],/NORMAL,COLOR=outsidecolor
        POLYFILL,!X.WINDOW[[0,1,1,0,0]],!Y.WINDOW[[0,0,1,1,0]],/NORMAL,COLOR=backgroundcolor
    ENDIF
ENDIF

IF N_ELEMENTS(normalonly) EQ 0 THEN normalonly=0
IF normalonly EQ 1 THEN noplot=1
IF N_ELEMENTS(noplot) EQ 0 THEN noplot=0
IF N_ELEMENTS(xtangle) EQ 0 THEN xtangle=0

IF N_ELEMENTS(leftaxis) EQ 0 THEN leftaxis=1
IF N_ELEMENTS(rightaxis) EQ 0 THEN rightaxis=1

IF noplot EQ 0 THEN BEGIN
    FOR i=0,n-1 DO BEGIN
; ----------- set up coordiante system -------
        xfrac=xxfrac[i:i+1]
        xframe=frame
        xframe([0,2])=xframe[0]+xfrac*width
        
        xrange=xxrange[i:i+1]*1.D0
        
        xnoerase=1
;        IF i EQ 0 THEN xnoerase=noerase

        plot,[0,1],/NODATA,NOERASE=xnoerase, $
          POSITION=xframe,/NORMAL, $
          YSTYLE=5, $
          XSTYLE=5, $
          XRANGE=xrange, $
          YLOG=ylog, $
          YRANGE=yrange


; --------- THE PLOT ------------
        IF N_ELEMENTS(proc) NE 0 THEN BEGIN
            CALL_PROCEDURE,proc,xrange,yrange
        ENDIF
        xmax=MAX(xrange,MIN=xmin) ; for decreasing x values
        IF oneline EQ 1 THEN BEGIN
            FOR ix=0,nx-1 DO BEGIN
                x1n=nn[ix]-1
                IF i EQ n-1 THEN BEGIN
                    end1[ix]=x1n
                ENDIF ELSE BEGIN            
                    IF xrange[0] LT xrange[1] THEN BEGIN
                        end1_temp=WHERE(x[0:x1n,ix] LE xmax)
                    ENDIF ELSE BEGIN
                        end1_temp=WHERE(x[0:x1n,ix] GE xmin)
                    ENDELSE
                    end1[ix]=end1_temp(N_ELEMENTS(end1_temp)-1)
                    IF end1[ix] EQ -1 THEN end1[ix]=last1[ix]
                ENDELSE
                IF last1[ix] LT end1[ix] THEN BEGIN
                    sel1=last1[ix]+1L+LINDGEN(MAX(end1[ix]-last1[ix]))
                    IF xmulti EQ 0 THEN BEGIN
                        iy0=0L
                        iy1=ny-1L
                    ENDIF ELSE BEGIN
                        iy0=ix
                        iy1=ix
                    ENDELSE
                    FOR iy=iy0,iy1 DO BEGIN
                        xyz=convert_coord(x[sel1,ix],y[sel1,iy] >0.5D0^(SGN(yrange[0]))*yrange[0] <2.D0^(SGN(yrange[1]))*yrange[1],/DATA,/TO_NORMAL,/DOUBLE)
                        y1[sel1,iy]=xyz[1,*]
                    ENDFOR
                    x1[sel1,ix]=xyz[0,*]
                ENDIF
                last1[ix]=end1[ix]
            ENDFOR
        ENDIF ELSE BEGIN
            IF nodata EQ 0 THEN BEGIN
                FOR iy=0,ny-1 DO BEGIN
                    IF xmulti EQ 1 THEN ix=iy ELSE ix=0L
                    IF (i EQ 0) THEN BEGIN
                        sel0=WHERE(x[0:nn[ix]-1,ix] GE xmin) 
                    ENDIF ELSE BEGIN
                        sel0=WHERE(x[0:nn[ix]-1,ix] GT xmin) 
                    ENDELSE
                    IF (i EQ n-1) THEN BEGIN
                        sel1=WHERE(x[sel0,ix] LE xmax) 
                    ENDIF ELSE BEGIN
                        sel1=WHERE(x[sel0,ix] LT xmax) 
                    ENDELSE
                    sel=sel0[sel1]
                    oplot,x[sel,ix],y[sel,iy], $
                      COLOR=color(iy), $
                      LINE=line(iy), $
                      THICK=thick(iy),NOCLIP=0
                ENDFOR
            ENDIF
        ENDELSE

        IF noframe EQ 0 THEN BEGIN
; ---------- x-axes ------------
            xdist=xrange(1)-xrange(0)
            xwidth=width*(xfrac(1)-xfrac(0))
            IF i LT n-1 THEN BEGIN
                xrange(1)=xrange(1)-xdist*delta/xwidth
                xframe(2)=xframe(2)-delta
            ENDIF
            IF i GT 0 THEN BEGIN
                xrange(0)=xrange(0)+xdist*delta/xwidth
                xframe(0)=xframe(0)+delta
            ENDIF

            xtickformat=xxtickformat[i]
            xlabel=''
            IF N_ELEMENTS(xxlabel) GE n THEN BEGIN
                xtickformat='spaghetti_nolabel'
                xlabel=xxlabel[i]
            ENDIF

            plot,[0,1],/NODATA, $
              POSITION=xframe,/NORMAL, $
              XSTYLE=1,YSTYLE=5, $
              XRANGE=xrange, $
              /NOERASE, $
              XTICKS=xxticks[i], $
              XTICKINTERVAL=xxtickinterval[i], $
              XTICKFORMAT=xtickformat, $
              XTITLE=xxtitle[i]

; --------- labels --------
            IF xtangle EQ -90 THEN BEGIN
                align=0
                xpos=frame[0]+0.5*(xxfrac[i]+xxfrac(i+1))*width-0.5D0*chsize[0]
                ypos=frame[1]-0.5*chsize[1]
                angle=-90
            ENDIF ELSE IF xtangle EQ 90 THEN BEGIN
                align=1
                xpos=frame[0]+0.5*(xxfrac[i]+xxfrac(i+1))*width+0.5D0*chsize[0]
                ypos=frame[1]-0.5*chsize[1]
                angle=90
            ENDIF ELSE BEGIN
                xpos=frame[0]+0.5*(xxfrac[i]+xxfrac(i+1))*width
                ypos=frame[1]-1.25*chsize[1]
                align=0.5
                angle=0
            ENDELSE
            XYOUTS,xpos,ypos,xlabel,ALIGNMENT=align,ORIENTATION=angle,/NORMAL
            
; --------- side axis setup --------
            plot,[0,1],/NODATA, $
              POSITION=frame,/NORMAL, $
              XSTYLE=5,YSTYLE=5, $
              XRANGE=xrange, $
              /NOERASE

; -------- left axis ---------
            IF (i EQ 0) AND (leftaxis NE 0) THEN BEGIN
                axis,YAXIS=0,YTITLE=ytitle, $
                  YLOG=ylog, $
                  YRANGE=yrange, $
                  YSTYLE=1
            ENDIF    

; -------- right axis ---------
            IF (i EQ n-1) AND (rightaxis NE 0)  THEN BEGIN
                axis,YAXIS=1, $
                  YLOG=ylog,YRANGE=yrange, $
                  YSTYLE=1, CHARSIZE=1.D-7
            ENDIF

            
; -------- interface ---------
            IF i LT n-1 THEN BEGIN
                xframe=frame
                xframe([0,2])=xframe(0)+xxfrac(i+1)*width

;        xpol=[-1,-1,0]*delta
;        ypol=[0,1,0]*bheight
                
;        polyfill,xpol+xframe(0),ypol+xframe(1),/NORMAL
;        plots,xpol+xframe(2),ypol+xframe(3),/NORMAL

; some eraseing
                xpol=[-1,-1,1]*delta
                ypol=[-1,1,-1]*bheight
                polyfill,xpol+xframe(0),ypol+xframe(1),/NORMAL,COLOR=outsidecolor

                xpol=[-1,1,1]*delta
                ypol=[1,-1,1]*bheight
                polyfill,xpol+xframe(2),ypol+xframe(3),/NORMAL,COLOR=outsidecolor

                xpol=[0, 1,1]*delta
                ypol=[0,-1,0]*bheight
                polyfill,xpol+xframe(0),ypol+xframe(1),/NORMAL,COLOR=xxedgecolor[i]

                xpol=[-1,-1,0]*delta
                ypol=[ 0, 1,0]*bheight
                polyfill,xpol+xframe(2),ypol+xframe(3),/NORMAL,COLOR=xxedgecolor[i]

; the "Z"        
                xpol=[-1,-1,1,1]*delta
                ypol=[0,1,-1,0]*bheight

                plots,xpol+xframe(0),ypol+xframe(1),/NORMAL
                plots,xpol+xframe(2),ypol+xframe(3),/NORMAL

; the tick mark on the break
                IF xxnotick[i] EQ 0 THEN BEGIN
                    xrange=xxrange(i+1)+[-1,1]
                    xframe=frame
                    xframe([0,2])=xframe(0)+xxfrac(i+1)*width+[-1,1]*1.D-7
                    
                    IF N_ELEMENTS(xxano) NE 0 THEN BEGIN
                        xano=xxano(i)
                    ENDIF ELSE BEGIN
                        v=xxrange(i+1)
                        s=STRTRIM(STRING(v),1)
                        s1=''
                        v1=0
                        i1=0
                        WHILE (ABS(v1-v) GT 1.D-7*ABS(v)) DO BEGIN
                            s1=s1+STRMID(s,i1,1)
                            v1=DOUBLE(s1)
                            i1=i1+1
                        ENDWHILE
                        xano=s1
                    ENDELSE
                    
                    plot,[0,1],/NODATA,/NOERASE, $
                      XRANGE=xrange, $
                      YRANGE=yrange, $
                      YSTYLE=4,XSTYLE=1,$
                      POSITION=xframe, $
                      xtickv=xxrange(i+1)+[-1,0,1]*1D-9, $
                      xticks=2, $
                      xtickname=[' ',xano,' '], $
                      xticklen=xxticklen[i]
                ENDIF        
            ENDIF
        ENDIF
    ENDFOR

    IF noframe EQ 0 THEN BEGIN
; --------- the x-title -------------
        xyouts,0.5D0*(frame(0)+frame(2)),frame0(1)+0.25*chsize(1),xtitle,ALIGN=0.5,/normal

; --------- the figure title -------------
        IF N_ELEMENTS(title) NE 0 THEN BEGIN
            xyouts,0.5D0*(frame[0]+frame[2]),frame0[3]-1.25*chsize[1],title,ALIGN=0.5,/normal
        ENDIF
    ENDIF
ENDIF

IF N_ELEMENTS(nonormal) EQ 0 THEN nonormal=0
IF nonormal EQ 1 THEN setnormal=0
IF N_ELEMENTS(setnormal) EQ 0 THEN setnormal=1
; --------- final setup for further plots --------
IF setnormal eq 1 THEN BEGIN
    plot,[0,1],/NODATA, $
      POSITION=frame,/NORMAL, $
      XSTYLE=5,YSTYLE=5, $
      XRANGE=xrange, $
      YRANGE=yrange, $
      YLOG=ylog, $
      /NOERASE
ENDIF

IF (oneline EQ 1) AND (noplot EQ 0) THEN BEGIN
    IF nodata EQ 0 THEN BEGIN
        FOR iy=0,ny-1 DO BEGIN
            IF xmulti EQ 1 THEN ix=iy ELSE ix=0
            x2=x1[0:nn[ix]-1,ix]
            y2=y1[0:nn[iy]-1,iy]
            IF reduce NE 0 THEN reduce,x2,y2,res=reduce
            plots,x2,y2, $
              NOCLIP=0, /NORMAL, $
              COLOR=color(iy), $
              LINE=line(iy), $
              THICK=thick(iy)
        ENDFOR
    ENDIF
    xmap=x1
    ymap=y1
ENDIF    

END

;-----------------------------------------------------------------------
;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
;-----------------------------------------------------------------------

PRO test

fac=1.D12

n=60000L
y=1.D0/(cos((2*!DPi*DINDGEN(n+1))/n)+1+1.D0/n)/n
x=DINDGEN(n+1)*fac


xxfrac=[0.0,0.3,0.5,0.7,1.0]
xxrange=[0,20000,29800,30200,60000D0]*fac
xxano=['a','b','c']
xxticklen=0

spaghetti,x,y, $
  XXFRAC=xxfrac, $
  XXRANGE=xxrange, $
  XXANO=xxano, $
  XXNOTICK=1, $
  YRANGE=[-0.1,1.1], $
  XXTICKINTERVAL=[5000,2000,50,5000]*fac, $
  XXTICKFORMAT="(I5)", $
  XXTITLE=['A','B','C','D'], $
    /ONELINE,XMAP=xmap,/NODATA,/NOFRAME,/NONORMAL
  
xyz=convert_coord(x,y,/DATA,/TO_NORMAL,/DOUBLE)

plots,xmap,xyz[1,*],/NORMAL,NOCLIP=0

spaghetti,x,y, $
  XXFRAC=xxfrac, $
  XXRANGE=xxrange, $
  XXANO=xxano, $
  XXNOTICK=1, $
  YRANGE=[-0.1,1.1], $
  XXTICKINTERVAL=[5000,2000,50,5000]*fac, $
;  XXTICKFORMAT="(I5)", $
  XXTITLE=['A','B','C','D'], $
    /NOERASE,/NODATA


END