% Look for "to make flexible" and "to change into" items, and "necessary?"
numeric temp_a[], temp_b;
numeric horizontal_span, vertical_span;
numeric no_of_instr, char_no, index, instr_, last_note[];
numeric vert_shift, instr_shift[];
numeric quanta[], width[], stretched_width[], to_stretch, last_stretch;
numeric note_axis[], ltx[], rtx[], rtx_[];          % rtx_ is the extreme of the current char.
numeric correction[];                               % x in the equations
numeric note, a[], b[], c[], d[], e[], f[], g[];
%numeric aggr_note[][];                           %   Highest (if #1 is 1) or lowest (#1=-1) note at instr_
numeric aggr_note[][][];
numeric instr_clef[], lowest_note[], highest_note[];
numeric stem_dir[][], stem_coef[], stem_extreme;
numeric i, j, s;                                   %   To change into i_, j_
numeric stem_axis[][];
numeric multi_beam_offset;
pair stem_point[][];
path temp_path;
picture note_[][];                              %   Old char_; [index][instr_]
picture pict_, _pict;
numeric next_level, curr_level;
string arg_text, temp_text;
string to_do, doing_, horizontal_elements[];
string curr_instr;
string staff_[];
numeric bar_type[];                             % Three-digit number. :|: is 111. | is (0)10.
numeric clef_sep, key_sep;
numeric total_natural_length;                   % In pt
numeric total_factors;                          % Total spacing factors (results of sf_, function of the quanta)
numeric spacing_unit;                           % What the stretching function multiplies. Default: a sharp's width
numeric least_spaced;                           % Quanta of the least-spaced rhythmic value. Usually 32 or 16
numeric no_of_blocks;                           % Number of blocks so far included in the line.
numeric index_offset;                           % To convert the first index of a line into 1
numeric no_of_beams;                            % Number of beams in the current line.
string beam_list[];                             % The different beams to be traced at the end of the line.
string tie_list[];
numeric no_of_ties, tie_[], open_tie[], closed_tie[];
string clef_changes;

%Initialization
to_do:="";
for octave=1 upto 8:
    g[octave]=3.5*(octave-4);
    a[octave]=g[octave]-3;
    b[octave]=g[octave]-2.5;
    c[octave]=g[octave]-2;
    d[octave]=g[octave]-1.5;
    e[octave]=g[octave]-1;
    f[octave]=g[octave]-.5;
endfor;
index:=0;
instr_:=0;
stem_length=3.5;
clef_sep=2interline;
key_sep=horizontal_axis;
horizontal_span:=10regular_width;
vertical_span:=10interline;
multi_beam_offset:=1/4;
note_axis0:=0;
width0 := 0;
stretched_width0:=0;
curr_level:=0;
next_level:=0;
total_natural_length:=0;
least_spaced:=32;
spacing_unit:=horizontal_axis+framingspace;
first_included:=1;
no_of_blocks := 0;
total_factors := 0;
next_char:=1;
no_of_beams := 0;
no_of_ties:=0;
index_offset := 0;

% Level-1 elements.

def add_halfheads(text nhs) =
    for s=nhs:
        note:=s+instr_clef[instr_];
        lowest_note[instr_] := min(lowest_note[instr_],note);
        highest_note[instr_] := max(highest_note[instr_],note);
        half_note_head shifted(0,interline*note);   %   To change into notehead.
%   Ledger lines follow. To make flexible.
        if (note>3.5) or (note<-1.5) :
            pickup penrazor xscaled line_thick rotated 90;
            for i:=(1+note/abs(note)*3) step (note/abs(note)) until note :
                draw (-.5regular_width,i*interline)--(.5regular_width,i*interline);
            endfor;
        fi;
    endfor;
    stem_axis[index][instr_] := xpart(stem_point[index][instr_]);
    find_extremes(nhs);                     % Vertical extremes
    enddef;

def add_noteheads(text nhs) =
    for s=nhs:
        note:=s+instr_clef[instr_];
        lowest_note[instr_] := min(lowest_note[instr_],note);
        highest_note[instr_] := max(highest_note[instr_],note);
        regular_note_head shifted(0,interline*note);   %   To change into notehead.
        note_axis[index] := 1/2regular_width;
        width[index] := 1/2regular_width;
%   Ledger lines follow. To make flexible.
        if (note>3.5) or (note<-1.5) :
            pickup penrazor xscaled line_thick rotated 90;
            for i:=(1+note/abs(note)*3) step (note/abs(note)) until note :
                draw (-.5regular_width,i*interline)--(.5regular_width,i*interline);
            endfor;
        fi;
    endfor;
    stem_axis[index][instr_] := xpart(stem_point[index][instr_]);
    find_extremes(nhs);                     % Vertical extremes
    enddef;
    
def add_rest(expr p, q) =
    addto currentpicture also rest_[p] shifted (0, q*interline);
    width[index] := rest_width[p];
    enddef;    

def add_wholeheads(text nhs) =
    for s=nhs:
        note:=s+instr_clef[instr_];
        lowest_note[instr_] := min(lowest_note[instr_],note);
        highest_note[instr_] := max(highest_note[instr_],note);
        whole_note shifted(0,interline*note);   %   To change into notehead.
%   Ledger lines follow. To make flexible.
        if (note>3.5) or (note<-1.5) :
            pickup penrazor xscaled line_thick rotated 90;
            for i:=(1+note/abs(note)*3) step (note/abs(note)) until note :
                draw (-.5regular_width,i*interline)--(.5regular_width,i*interline);
            endfor;
        fi;
    endfor;
    stem_axis[index][instr_] := xpart(stem_point[index][instr_]);
    find_extremes(nhs);                     % Vertical extremes
    enddef;


def barline_(expr number) =
    %The first digit indicates pre-barline stuff.
    temp:=floor(number/100);
    %The middle digit indicates the barline itself.
    temp:=floor(number/10 - 10temp);
    if temp = 1 :
        pickup penrazor xscaled line_thick; 
        draw (0, min(scantokens staff_[instr_])*interline)--(0, max(scantokens staff_[instr_])*interline);
%        rtx_[instr_]:=-1/2regular_width;
        stem_dir[index][instr_] := 0;
    elseif temp=2 : 
        pickup penrazor xscaled line_thick;
        draw (0, min(scantokens staff_[instr_])*interline)--
            (0, max(scantokens staff_[instr_])*interline);
        pickup penrazor xscaled 1/2interline;
        draw (1/2interline+1/2line_thick, min(scantokens staff_[instr_])*interline)--
            (1/2interline+1/2line_thick, max(scantokens staff_[instr_])*interline);
%        rtx_[instr_]:=-1/2regular_width;
%        ltx[instr_]:=interline-1/2regular_width;
        width[index] := 3/4interline + 1/2line_thick;
        stem_dir[index][instr_] := 0;
    fi;
    %The last digit indicates post-barline stuff.
    enddef;

% Level-2 elements.

def add_dot (suffix nhs) = 
    temp_text:= str nhs;
    if curr_level < 2 :                                %   Deferring...
        to_do := if next_level>0 : to_do & fi curr_instr & "add_dot(" & temp_text & ");";  
        curr_instr := "";
        next_level:=2;
    else:                                           %   Doing...
        for i:=nhs :
            note:=floor(i+instr_clef[instr_])+.5;
            pickup pencircle scaled 1/3interline;
            drawdot (1/2regular_width+framingspace, note*interline);
        endfor;
    fi;
    enddef;

def add_flag (suffix number) =            
    temp_text := str number;
    if curr_level < 2 :                                %   Deferring...
        to_do := if next_level>0 : to_do & fi curr_instr & "add_flag(" & temp_text & ");";  
        curr_instr := "";
        next_level:=2;
    else:                                           %   Doing...
    temp_a := scantokens temp_text;
    if number = 1 :
        eighth_flag((0, interline*aggr_note[index][instr_][-stem_coef[instr_]]))
            if stem_coef[instr_]=1 : 
                reflectedabout ((stem_axis[index][instr_],0),(stem_axis[index][instr_],1)); 
                rtx_[instr_]:=max(rtx_[instr_],horizontal_axis);
            else:
                rotatedaround ((0, interline*aggr_note[index][instr_][-stem_coef[instr_]]),180)
            fi;
    fi;
    fi;
    enddef;

def regular_stem =
    if curr_level<2:                                %   Deferring...
        to_do := if next_level>0 : to_do & fi curr_instr & "regular_stem;";  
        curr_instr := "";
        next_level:=2;
    else:                                           %   Doing...
    if unknown stem_dir[index][instr_] : 
        stem_dir[index][instr_] = stem_coef[instr_];
    fi;
        pickup penrazor xscaled line_thick;
        draw (stem_axis[index][instr_],interline*aggr_note[index][instr_][stem_dir[index][instr_]])
            --(stem_axis[index][instr_], 
                interline*(aggr_note[index][instr_][stem_dir[index][instr_]]+stem_dir[index][instr_]*stem_length));
    fi; 
    enddef;
    
% Level-3 elements.

def add_flat(suffix notext) =
    temp_text:=str notext;
    if curr_level<3 :
        to_do := if next_level>0 : to_do & fi curr_instr & "add_flat(" & temp_text & ");";
        curr_instr := "";
        next_level :=3;
    else:
        temp_a:= scantokens temp_text;
        addto currentpicture also flat shifted (-1/2regular_width-ltx[instr_],interline*(_note(temp_a,instr_,index)));
        ltx[instr_] := ltx[instr_]+horizontal_axis;
    fi;
    enddef;
    
def add_natural(suffix notext) =
    temp_text:=str notext;
    if curr_level<3:                                %   Deferring...
        to_do := if next_level>0 : to_do & fi curr_instr & "add_natural(" & temp_text & ");";  
        curr_instr := "";
        next_level:=3;
    else:                                           %   Doing...
        temp_a:= scantokens temp_text;        
        addto currentpicture also natural shifted (-1/2regular_width-ltx[instr_],interline*(_note(temp_a,instr_,index)));
        ltx[instr_] := ltx[instr_]+horizontal_axis;
    fi;
    enddef;

def add_sharp(suffix notext) =
    temp_text:=str notext;
    if curr_level<3:                                %   Deferring...
        to_do := if next_level>0 : to_do & fi curr_instr & "add_sharp(" & temp_text & ");";  
        curr_instr := "";
        next_level:=3;
    else:                                           %   Doing...
        temp_a:= scantokens temp_text;        
        addto currentpicture also sharp shifted (-1/2regular_width-ltx[instr_],interline*(_note(temp_a,instr_,index)));
        ltx[instr_] := ltx[instr_]+horizontal_axis-framingspace;
    fi;
    enddef;

% Level-4 elements.

% Level-5 elements.

def beam(text t_) = 
    no_of_beams := no_of_beams + 1;
    beam_list[no_of_beams] := decimal(instr_); 
    for i:= t_ :
        beam_list[no_of_beams] := beam_list[no_of_beams] & "," & decimal(i);
    endfor;
    enddef;

def beams_ = 
    for k := 1 upto no_of_beams :    
        temp_text:="extract_beam_data(" & beam_list[k] & ");"; 
        scantokens temp_text;
        beam_list[k] := "";
        prepare_beam;
        currentpicture:=note_[closest_index][instr_];
        draw_beam;
        draw_stems;
        note_[closest_index][instr_]:=currentpicture;
        currentpicture:=note_[index][instr_];
    endfor;
    enddef;

def add_tie (text nhs) = 
    forsuffixes s:=nhs :
        tie_list[instr_] := if known tie_list[instr_] : tie_list[instr_] & "," & fi decimal index & "," & str s;
    endfor;
    enddef;

def close_tie (text nhs) =
    if index  > 1 :
        forsuffixes s:=nhs :
            tie_list[instr_] :=  tie_list[instr_] & ",-" & decimal index & "," & str s;
        endfor;
    fi;
    enddef;

def ties_ =  
    for i := 1 upto no_of_instr :
        if known tie_list[i] : 
            k:=100;
            for j:= scantokens tie_list[i] :
                if k=100 :        % j is an index
                    if j>0 :    % this is an opening index
                        no_of_ties:=no_of_ties+1;
                        open_tie[no_of_ties]:=j;
                        k:=200;
                    else:       % this is a closing index
                        k:=-j;
                    fi;
                elseif k=200 :        % j is an opening note
                    tie_[no_of_ties] := j;
                    k:=100;
                else:           % j is a closing note
                    for l:=1 upto no_of_ties :
                        if tie_[l] = j :    % this is it
                            closed_tie[l] := k;
                        fi;
                    endfor;
                    k:=100;
                fi;
            endfor;
            for j:=1 upto no_of_ties :
                if known closed_tie[j] :
                    draw_tie(i, open_tie[j],closed_tie[j],_note(tie_[j],i,open_tie[j]));
                else:
                    draw_open_tie(i, open_tie[j],_note(tie_[j],i,open_tie[j]));
                fi;
            endfor;
            numeric tie_[], open_tie[], closed_tie[];
            no_of_ties:=0;
        fi;
    endfor;
    enddef;

def draw_tie(expr ins_, oi_, ci_, tn_) =
    tie_coef := -stem_dir[oi_][ins_];
    begingroup; save x, y;
    path p;
    clearit;
    z0 = (0, interline*(tn_));
    x2=total_width(oi_) 
        for m:=oi_+1 upto ci_-1 : + note_axis[m] + total_width(m) endfor 
        + note_axis[ci_];
    y2=y0;
    y1=y0+tie_coef*interline; 
    x1 = .5[x0,x2];
    p = z0..{right}z1..z2;
    z3 = p intersectionpoint 
        ((1/2horizontal_axis, interline*(tn_-1))--(1/2horizontal_axis, interline*(tn_+1)));
    z4 = p intersectionpoint
        ((x2-1/2horizontal_axis,interline*(tn_-1))--(x2-1/2horizontal_axis,interline*(tn_+1)));
    penpos3(line_thick, 90+notehead_angle);
    penpos1(2line_thick, 90);
    penpos4(line_thick, 90-notehead_angle);
    penstroke z3e..{right}z1e..z4e;
    addto note_[oi_][ins_] also currentpicture 
        shifted (tie_coef*.5framingspace, tie_coef*2line_thick);%*(tn_-floor(tn_)));
    clearit;
    endgroup;
    enddef;

def draw_open_tie(expr ins_, oi_, tn_) =
    tie_coef := -stem_dir[oi_][ins_];
    begingroup; save x, y;
    path p;
    clearit;
    z0 = (0, interline*(tn_));
    x2=total_width(oi_) 
        for m:=oi_+1 upto index : + note_axis[m] + total_width(m) endfor
        if (tn_ < max(scantokens staff_[instr_]) + 1) and 
            (tn_ > min(scantokens staff_[instr_]) - 1) :
            + regular_width 
        fi;
    y2=y0;
    y1=y0+tie_coef*interline;
    x1 = .5[x0,x2];
    p = z0..{right}z1..z2;
    z3 = p intersectionpoint 
        ((1/2horizontal_axis, interline*(tn_-1))--(1/2horizontal_axis, interline*(tn_+1)));
    z4 = p intersectionpoint
        ((x2-1/2horizontal_axis,interline*(tn_-1))--(x2-1/2horizontal_axis,interline*(tn_+1)));
    penpos3(line_thick, 90+notehead_angle);
    penpos1(2line_thick, 90);
    penpos4(line_thick, 90-notehead_angle);
    penstroke z3e..{right}z1e..z4e;
    addto note_[oi_][ins_] also currentpicture 
        shifted (framingspace, tie_coef*2line_thick);%*(tn_-floor(tn_)));
    clearit;
    endgroup;
    enddef;

% Functions

def clef_ (expr number) =
    clef_changes := decimal(index) & "," & decimal(instr_) & "," & decimal(instr_clef[instr_])
        if known clef_changes : & "," & clef_changes fi;
    instr_clef[instr_]:=number;
    small_clef(number);
    enddef;

def small_clef (expr number) = 
    if curr_level < 4 :                             %  Deferring...
        to_do := if next_level>0 : to_do & fi curr_instr & "small_clef(" & decimal(number) & ");"; 
        curr_instr := "";
        next_level:=4;
    else:                                           %   Doing...
        addto currentpicture also
            if instr_clef[instr_] = 0 :
                small_treble_clef shifted (-2.8interline-1/2regular_width-ltx[instr_],0);
                ltx[instr_] := ltx[instr_] + 2.4interline+1/2regular_width + framingspace;
            elseif instr_clef[instr_] = 6 :
                small_bass_clef shifted (-2.75interline-1/2regular_width-ltx[instr_],0);
                ltx[instr_] := ltx[instr_] + 2.4interline + 1/2regular_width + framingspace;
%                pickup pencircle scaled line_thick;
%                draw (-ltx[instr_],-interline)--(-ltx[instr_],3interline);
            fi;
    fi;
    enddef;    

def draw_beam =
    hor_start := stem_axis[index][instr_]+.5line_thick 
        for i:=closest_index upto (index-1): + note_axis[i+1] + total_width(i) endfor;
    pickup penrazor xscaled .5interline rotated 90;
    for i:=beam_elements downto 2 :
        hor_end := hor_start 
            for j:=beam_index[i-1] upto (beam_index[i]-1) : - total_width(j) - note_axis[j+1] endfor 
            -stem_axis[index][instr_] + stem_axis[beam_index[i]][instr_] -line_thick;
        for j := 1 upto beam_no[i] :
            draw ((hor_start,interline*(beam_lift - beam_coef*(3j-2)/4))--
                (hor_end,interline*(beam_lift - beam_coef*(3j-2)/4))) transformed yslanted;
        endfor;
% An extra parameter for each note of the stem will indicate whether it's + or - horizontal_axis
        if beam_no[i] > beam_no[i+1] :
            for j :=2 upto beam_no[i] :
                draw ((hor_start,interline*(beam_lift - beam_coef*(3j-2)/4))--
                    (hor_start-horizontal_axis,interline*(beam_lift - beam_coef*(3j-2)/4))) transformed yslanted;
            endfor
        fi;
        for j:= (beam_no[i-1]+1) upto beam_no[i] :
            undraw ((hor_start-line_thick,interline*(beam_lift - beam_coef*(3j-2)/4))--
                (hor_end+line_thick,interline*(beam_lift - beam_coef*(3j-2)/4))) transformed yslanted;
        endfor;
        highest_note[instr_]:=max(highest_note[instr_], ypart((
            if beam_dir < 0 : hor_end else: hor_start fi
            , beam_lift*interline) transformed yslanted)/interline);
        lowest_note[instr_]:=min(lowest_note[instr_], ypart((
            if beam_dir < 0 : hor_end else: hor_start fi
            , beam_lift*interline) transformed yslanted)/interline);
        hor_start := hor_end + line_thick;
    endfor;
    enddef;

def draw_stems =
    pickup penrazor xscaled line_thick;
    hor_start := stem_axis[index][instr_] 
        for i:=closest_index upto (index-1): + note_axis[i+1] + total_width(i) endfor;
    for i:=beam_elements downto 1 :
        draw (xpart(stem_point[beam_index[i]][instr_])
            for j:=closest_index upto (beam_index[i]-1): + note_axis[j+1] + total_width(j) endfor
            for j:=beam_index[i] upto (closest_index-1): - note_axis[j+1] - total_width(j) endfor,
            ypart(stem_point[beam_index[i]][instr_])
            + interline * if beam_coef =1 : min else: max fi (beam_aggr[i]1, 
                for j:=1 upto beam_aggr[i]0 : , beam_aggr[i][j] endfor))
            --((xpart(stem_point[beam_index[i]][instr_])
            for j:=closest_index upto (beam_index[i]-1): + note_axis[j+1] + total_width(j) endfor
            for j:=beam_index[i] upto (closest_index-1): - note_axis[j+1] - total_width(j) endfor,
            interline*beam_lift) transformed yslanted);
    endfor;
    enddef;

def end_of_block =
    temp:=quanta[index];
    quanta[index]:=0;
    if unknown desired_width :              % A \musicbox
        desired_width:=ideal_width;
    else :
        no_of_blocks := no_of_blocks + 1;
        if abs(desired_width-ideal_width) > abs(desired_width - ideal_width*(no_of_blocks+1)/no_of_blocks) :      % Don't make the line.
        else :
            make_line(next_char);
        fi
    fi;
    quanta[index]:=temp;
    enddef;
    
def ideal_width =
    (total_natural_length + (
    for i := 1 upto index-1 :
        max(sf_(quanta[i])*spacing_unit/pt-correction[i+1]/pt,0) +
    endfor
    sf_(quanta[index])*spacing_unit/pt))
    enddef;

def sf_(expr number) =
%    if number 1.618 ** (mlog(if number = 0 : 1 else : number/least_spaced fi)/mlog(2))
    if number = 0 : 0
    else :
        1.618 ** (mlog(abs(number)/least_spaced)/mlog(2))
    fi
    enddef;

def stretch_chars =
    for i := 1 upto index :
        total_factors := total_factors + sf_(quanta[i]);
    endfor;
    forever: 
        temp := 0;
        note_axis[index+1]:=0;
        correction[index+1]:=0;
        for i := 1 upto index :
            stretched_width[i] := 
                max(sf_(quanta[i])*(desired_width-total_natural_length)/total_factors
                    - correction[i+1]/pt,0);
            stretched_width[i] := stretched_width[i]*pt;
            temp := temp + (note_axis[i] + total_width(i))/pt;
        endfor;
        exitif (desired_width - temp >= 0) and (desired_width - temp <=0.005);
        total_natural_length := total_natural_length - (desired_width - temp);
    endfor;
    enddef;
    
def end_of_char =
    note_[index][instr_]:=currentpicture;
    clearit;
    ensure_info;
    curr_level := next_level;
    next_level := 0;
    doing_ := to_do & "end_of_char;";
    to_do := "";
    instr_:=0;
    if curr_level = 0 :                                                 %   Which means last next_level=0
        finalize_char; 
    else:
        scantokens doing_;
    fi;
    enddef;

def ensure_info =
    if ((next_level=0) and (curr_level<6)) or (next_level=6) : 
        if index = 1 :
            for i:= 1 upto no_of_instr :
                note_axis1 := max(note_axis1, ltx[i]);
            endfor;
        else:
            for i:=1 upto no_of_instr :
                if totalweight note_[index][i]>0 :
                    avlbl_space := 0 for j:=last_note[i]+1 upto index-1 : 
                        + note_axis[j] + width[j] endfor + note_axis[index];
                    if avlbl_space < (rtx[i] + ltx[i]) : 
                        find_new_axis(i); 
                        fi;
            fi;
            endfor;
        fi;
    fi;
    enddef;

def extract_beam_data(text t) =
    j:=0; 
    numeric beam_elements, beam_index[], beam_no[], beam_aggr[][];
    beam_index0:=0;
    beam_elements:=0;
    instr_ := 0;
    for i:=t :
        if instr_ = 0 : instr_:=i;
        else :
            j:=j+1;
            if j=1 :        % New element, index
                beam_elements:=beam_elements+1;
                beam_index[beam_elements]:=i-index_offset;
            elseif j=2 :    % No. of beams
                beam_no[beam_elements]:=i;
            elseif j<1 :    % One of the notes in the aggr. Position 0 in the array is no. of notes.
                beam_aggr[beam_elements][1-j] := _note(i, instr_, beam_index[beam_elements]);
            else :          % Number of notes
                beam_aggr[beam_elements][0]:=i;
                j:=-i;
            fi;
        fi;
    endfor;
    beam_no[beam_elements+1]:=0;
    enddef;
    
def finalize_char =
    for i:=1 upto no_of_instr :
        if totalweight(note_[index][i])>0 : last_note[i]:=index; rtx[i]:=rtx_[i]; fi;
%        prev_staff(i);
    endfor;
    find_pos;
    total_natural_length:=total_natural_length+(note_axis[index]+width[index])/pt;
    curr_level:=0;
    enddef;

def find_extr(expr p_, coeff_) =
    clearit;
    tx:= coeff_ * horizontal_span;
    fx:=0;
    pickup pencircle scaled epsilon
    forever: xtx:=.5[tx,fx]; exitif abs(tx-fx)<=.1;
        currentpicture:=p_;
        undraw(xtx,-vertical_span)--(xtx,vertical_span);
        if totalweight(currentpicture)>=totalweight(p_) : tx else: fx fi 
        :=xtx;
    endfor;
    xtx[instr_]:= coeff_*xtx;
    clearit;
    enddef;
    
def find_extremes(text notes_) =              %   Will be a text expresion, this macro decides everything
    aggr_note[index][instr_][1]:=-100;             % Top note
    aggr_note[index][instr_][-1]:=100;            % Bottom note
    for i:=notes_ :
        aggr_note[index][instr_][1]:=max(aggr_note[index][instr_][1],i+instr_clef[instr_]);
        aggr_note[index][instr_][-1]:=min(aggr_note[index][instr_][-1],i+instr_clef[instr_]);
    endfor;
    stem_coef[instr_]:= if (aggr_note[index][instr_][1]+aggr_note[index][instr_][-1])>=2: - fi 1;
    enddef;

def find_inclination (text options) =
    yslanted := identity;
    beam_height := 0;
    pair beam_pt[];     % The definition of beam_pt's has to happen even if there's no decision to be made,
                        % because find_position uses beam_pt1.
    beam_pt[closest_no] := (0,0);
    for i:=closest_no-1 downto 1:
        xpart beam_pt[i] = xpart(beam_pt[i+1])-stem_axis[beam_index[i]][instr_]-note_axis[beam_index[i]]
            for j:=beam_index[i] upto (beam_index[i+1]-1) :
                - note_axis[j] - total_width(j) endfor
            + note_axis[beam_index[i]] + stem_axis[beam_index[i]][instr_];
        ypart beam_pt[i] = interline * (beam_aggr[i]1-closest_note);
    endfor;
    for i:=closest_no+1 upto beam_elements:
        xpart beam_pt[i] = xpart(beam_pt[i-1]) - stem_axis[beam_index[i]][instr_] - note_axis[beam_index[i]]
            for j:=beam_index[i-1] upto (beam_index[i]-1) :
                + note_axis[j] + total_width(j) endfor
            + note_axis[beam_index[i]] + stem_axis[beam_index[i]][instr_];
        ypart beam_pt[i] = interline * (beam_aggr[i]1-closest_note);
    endfor;
    if beam_dir <> 0 :
        transform yslanted[];
        total_beam_width:=xpart(beam_pt[beam_elements])-xpart(beam_pt1);
        temp_a := 0;
        for i:= 1 upto beam_elements :
            temp_a := temp_a - beam_coef*ypart(beam_pt[i]);
        endfor;
%if beam_index[1]=1 :
%    show beam_dir;
%    for l:=1 upto beam_elements: 
%        show beam_pt[l];
%    endfor;
%fi;
        for i:= scantokens options :
            temp_b := -i*xpart(beam_pt1)/4total_beam_width; 
            (0,1) transformed yslanted[i] = (0,1);
            (xpart(beam_pt1),0) transformed yslanted[i] = (xpart(beam_pt1),interline*temp_b);
            (xpart(beam_pt[beam_elements]),0) transformed yslanted[i] = 
                (xpart(beam_pt[beam_elements]),interline*(temp_b-i/4));
            temp_a[i]:=0;
            for j := 1 upto beam_elements :
                if ypart (-beam_coef*(beam_pt[j] transformed yslanted[i])) >= -0.001 :
                    temp_a[i] := temp_a[i] - beam_coef*ypart(beam_pt[j] transformed yslanted[i]);
                else: 
                    temp_a[i] := 1000;
                fi;
                exitif temp_a[i] = 1000;
            endfor;
            if temp_a[i] - temp_a < .05 :
                transform yslanted;
                (0,1) transformed yslanted = (0,1);
                (-1,0) transformed yslanted = (-1,0) transformed yslanted[i] reflectedabout((0,0),(1,0));
                (1,0) transformed yslanted = (1,0) transformed yslanted[i] reflectedabout((0,0),(1,0));
                temp_a := temp_a[i];
                beam_height := -i/4; 
            fi;
        endfor;
    fi;
    enddef;
        
def find_ltx(expr p_) = find_extr(p_,-1); enddef;

def find_new_axis (expr inst) =
show index;
    pict_:=note_[last_note[inst]][inst] shifted %
        (width[last_note[inst]] for i:=(last_note[inst]+1) upto (index-1):
            -note_axis[i]-width[i] endfor,0);
    _pict:=note_[index][inst];
    cull _pict keeping (1,infinity);
    cull pict_ keeping (1,infinity);
    tx:=rtx[inst]+ltx[inst]+note_axis[index];
    fx:=ltx[inst];
%fx will be customized so that the user can indicate what fraction of ltx[inst] has to be taken into account.
%In any case, fx (the minimum new axis), being exactly ltx, makes room for the whole width of the previous note.
%Which is what is needed, because the stems are not yet drawn.
    nax:=tx;
    forever:
        clearit;
        exitif abs(tx-fx)<=.01;
        currentpicture := pict_;
        addto currentpicture also _pict shifted (nax,0);
        cullit;
        if totalweight(currentpicture)<totalweight(pict_)+totalweight(_pict):
            fx else: tx fi :=nax;
        nax:=.5[tx,fx]; 
    endfor;
    correction[index] := max(note_axis[index],nax) - note_axis[index];
    note_axis[index] := note_axis[index] + correction[index];
    enddef;

def find_pos =
%    last_stretch:=sfcode
%        *to_stretch*regular_width/total_space_factors;
    stretched_width[index]:=0;
    enddef;

def find_position =
    temp_a := closest_note + beam_coef*(stem_length-1/2);
    (ypart(beam_pt[beam_elements] transformed yslanted)-ypart(beam_pt[beam_elements]))/interline + beam_lift = 
        round(temp_a+.01beam_coef) + beam_coef/4;
    enddef;
    
def instr_no(expr number) =
    no_of_instr:=number;
    for i:=1 step 1 until number:
        highest_note[i]:=3;
        lowest_note[i]:=-1;
        last_note[i]:=0;
        rtx[i]:=0;
        note_0[i]:=nullpicture;
        instr_clef[i]:=0;
    endfor; 
    enddef;

def key_sign = relax
    if curr_level < 4 :
        to_do := if next_level>0 : to_do & fi "key_sign;"; 
        curr_instr := "";  % Maybe unnecessary
        next_level := 4;
    else:
    % Key signatures
        note_[index][instr_] := currentpicture;
        temp_text := "7";
        if key_/abs(key_+.1) = 1 :                     % Sharps
            relax;
        else :                                      % Flats
            j := 1; 
            for i := 1 upto -key_ :
                temp_text := decimal j & "," & temp_text;
                j :=  if odd(i): 1.5 else: -2 fi + j;
            endfor;
            for i := 1 upto no_of_instr :
                ltx[i] := ltx[i] +key_sep;
                for s:= scantokens(temp_text) :
                    exitif s = 7;
                    addto note_[index][i] also flat
                        shifted (-ltx[i]-horizontal_axis, 
                            interline * (s if instr_clef[i]=6:-1fi));
                    ltx[i] := ltx[i] + horizontal_axis;
                endfor;
            endfor;
        fi;
    fi;
    enddef;

def normal_clefs =
    if curr_level < 4 :                             %   Deferring...
        to_do := if next_level>0 : to_do & fi "normal_clefs;";  
        curr_instr := "";    % Maybe unnecessary
        next_level:=4;
    else:                                           %   Doing...
    %Clefs. To do: the font has to draw the clefs from 'trebleaxis' and 'bassaxis'
        note_[index][instr_] := currentpicture;
        for i := 1 upto no_of_instr :
            addto note_[index][i] also
                if instr_clef[i] = 0 :
                    treble_clef shifted (-3interline - clef_sep - ltx[i],0);
%                    note_axis[index] := max(note_axis[index], 3interline + clef_sep);
                    ltx[i] := ltx[i] + 3interline + clef_sep;
                elseif instr_clef[i] = 6 :
                    bass_clef shifted (-2.7interline - clef_sep - ltx[i],0);
%                    note_axis[index] := max(note_axis[index], 2.7interline + clef_sep);
                    ltx[i] := ltx[i] + 2.7interline + clef_sep;
                fi;
        endfor;
    fi;
    enddef;
    
def make_line(expr number) =
    if index > 0:
        quanta[index]:=0;
        stretch_chars;
        instr_shift1:=0;
        begingroup; save index;
        beams_;
        ties_;
        endgroup;
        for i:=2 upto no_of_instr:
            instr_shift[i]:=4+instr_shift[i-1]+highest_note[i-1]-lowest_note[i];
        endfor;
        line_height:=(instr_shift[no_of_instr]+highest_note[no_of_instr]+2)*interline#;
        line_depth:=-(lowest_note1-1)*interline#;
        for i:= 1 upto index:
            for j:=1 upto no_of_instr :
                clearit;
                staff_lines(i,scantokens staff_[j]);
                addto note_[i][j] also currentpicture;
            endfor;
            beginchar(number+i-1,(note_axis[i]+width[i]+stretched_width[i])*pt#/pt,line_height,line_depth);
                clearit;
            if no_of_instr > 1 :
                if i=1 :        %   Bracket
                    pickup penrazor xscaled line_thick;
                    draw (-1/2line_thick,-interline)--(-1/2line_thick,interline*(instr_shift2+3));
                    pickup pencircle xscaled 1/2interline yscaled line_thick rotated (-90-notehead_angle);
                    top rt z0 = (-2framingspace, interline*(instr_shift2+3));
                    bot lft z1 = (-horizontal_axis, interline*(instr_shift2/2+1)-1/2line_thick);
                    z2 = .5[z0,z1];
                    draw z0{dir(-90-notehead_angle)}..z2..{dir(-135)}z1;
                    pickup pencircle xscaled 1/2interline yscaled line_thick rotated (90+notehead_angle);
                    draw (z0{dir(-90-notehead_angle)}..z2..{dir(-135)}z1) 
                        reflectedabout ((0,interline*(instr_shift2/2+1)),(1,interline*(instr_shift2/2+1)));
                fi;
                if quanta[i]<=0 :            %   Inter-staff barline
                    pickup penrazor xscaled line_thick;
                    draw (note_axis[i],3interline)--(note_axis[i],interline*(instr_shift2-1));
                fi;
            fi;
                for j:=1 upto no_of_instr:
                    addto currentpicture also
                        note_[i][j] shifted (note_axis[i],interline*instr_shift[j]-(-1**i)*0/2line_thick);
                    note_[i][j]:=nullpicture;
                endfor;
            endchar;
        endfor;
    fi;
    next_char:= index + number;
    index_offset := index_offset+index;
    index:=0; 
    total_natural_length:=0;
    total_factors := 0;
    no_of_blocks := 0;
    no_of_beams := 0;
    string beam_list[], tie_list[];
    numeric stem_dir[][];
    string clef_changes;
    clearit;
showstats;
    enddef;

def time_sign = relax enddef;

def new_char(expr number) =
    curr_level:=1;
    next_level:=0;
    if index = 0 : 
        to_do := to_do & "time_sign;" & "key_sign;" & "normal_clefs;";
        next_level := 4;
    fi; 
    index:=index+1; 
    quanta[index]:= if number = 0 : least_spaced else : number fi;
%    total_factors := total_factors + abs(sf_(number));        
    if (number < 17) and (number > 0) : least_spaced:=16 fi;
    width[index]:=0;
    note_axis[index]:= 0;
    correction[index] := 0;
    for i:=1 upto no_of_instr:
        rtx_[i] := 0;
        ltx[i]:=0;
        note_[index][i]:=nullpicture;
    endfor;
    enddef;

def _note(expr n_, i_, _ix) =            % Note, instrument, index
    hide (
    _note_:=n_+instr_clef[i_];
    if known clef_changes :
        k:=0;
        for temp:=scantokens clef_changes :
            if k=0 :                    % An index
                k:= if temp<=_ix : - fi 1;
            elseif k=1 :                % An instrument
                if temp<i_ : 
                    k:=2;
                elseif temp=i_ :
                    k:=3;
                else:
                    k:=4;
                fi;
            elseif k=2 :                % Good index, but wrong instr
                k:=0;
            elseif k=3 :                % This is it
                _note_:=n_+temp;
                k:=0;
            else :                      % Good index, missed instr
                k:=0;
            fi;
            exitif k=-1;
        endfor;
    fi;
    )
    _note_
    enddef;
    

def prepare_beam =
% Direction of stems
    extr_beam_note1 := -100;
    extr_beam_note[-1] := 100;
    for i:=1 upto beam_elements :
        for j:=1 upto beam_aggr[i][0] :
            extr_beam_note1:=max(extr_beam_note1,beam_aggr[i][j]); % Top note
            extr_beam_note[-1]:=min(extr_beam_note[-1],beam_aggr[i][j]); % Bottom note
        endfor; 
    endfor; 
    beam_coef := if (extr_beam_note1-1) >= (1-extr_beam_note[-1]) : - fi 1; % Direction of stems
    for i:=1 upto beam_elements :
        if unknown stem_dir[beam_index[i]][instr_] : 
            stem_dir[beam_index[i]][instr_] = beam_coef;
        fi;
    index := beam_index[beam_elements];
    endfor;
% General direction (depending on the farthest-away note-heads in the first and last notes)
    temp_a := beam_aggr[1][1];
    temp_b := beam_aggr[beam_elements]1;
    for i:=2 upto beam_aggr[1][0] :
        temp_a := if beam_coef=1 : min else: max fi (temp_a,beam_aggr1[i]);
    endfor;
    for i:=2 upto beam_aggr[beam_elements]0 :
        temp_b := if beam_coef=1 : min else: max fi (temp_b,beam_aggr1[i]);
    endfor; 
    beam_dir := if temp_a = temp_b : 0* elseif temp_a > temp_b : - fi 1;
% Reduction to relevant notes (irrelevant notes are put in positions 2 and up of beam_aggr). 
% and finding of the closest note, having any inner beams into account too.
    closest_note:=-100beam_coef;
    for i:=if (beam_coef*beam_dir) >= 0 : 1 upto beam_elements else: beam_elements downto 1 fi:
        if beam_aggr[i]0 > 1 :
            if (beam_coef*beam_aggr[i]1) < (beam_coef*beam_aggr[i]2) :
                temp_a := beam_aggr[i]1; 
                beam_aggr[i]1 := beam_aggr[i]2;
                beam_aggr[i]2 := temp_a;
            fi;
        fi;
        temp_b := if beam_coef = 1 : max else: min fi 
            (closest_note,beam_aggr[i]1+beam_coef*beam_no[i]/100);
        if closest_note <> temp_b : 
            closest_index := beam_index[i]; closest_no:=i; fi;
        closest_note := temp_b;
    endfor;
    closest_note:=round(10temp_b)/10;
% "Mode" (careful or not)
    numeric beam_lift;
    temp_a := closest_note + beam_coef * (stem_length - 3*(beam_no[closest_no]-1)/4);
    transform yslanted;
    if (temp_a < max(scantokens staff_[instr_]) + 1) and (temp_a > min(scantokens staff_[instr_]) - 1) :
        "careful beam";
        find_inclination("beam_dir");
        find_position;
    else :
        find_inclination("2beam_dir, beam_dir");
        beam_lift := closest_note + beam_coef * (1/4 + stem_length);
    fi;        
    enddef;
        
def prev_staff (expr number) =
    clearit;
    staff_lines(index-1,scantokens staff_[number]);
    addto note_[index-1][number] also currentpicture;
    clearit;
    enddef;

def staff_lines(expr char_ind)(text t) = 
    pickup pencircle scaled line_thick;
    for j=t:
        draw ((0,j*interline)--(note_axis[char_ind]+width[char_ind]+stretched_width[char_ind],j*interline)) %
            shifted (-note_axis[char_ind],0);
    endfor;
    enddef;

def this_staff =
    for i:=1 upto no_of_instr:
        clearit;
        staff_lines(index,scantokens staff_[i]);
        addto note_[index][i] also currentpicture;
        clearit;
    endfor;
    enddef;

def total_width(expr idx_) =
    (width[idx_] + stretched_width[idx_])
    enddef;