CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C contour - melodic accent of a sequence of pitches
C   by Andy Allinger, 2015, released to the public domain
C   This program may be used by any person for any purpose.
C
C refer to:
C     "Melodic accent:  Experiments and a tentative model"
C     Joseph M. Thomassen
C     Journal of the Acoustical Society of America
C     v.71 n.6 pp.1596-1605, June 1982
C
C     There is an erratum for this article, which does not seem
C to be pertinent to the present computation.
C
C                          I/O LIST
C__Name______Type________I/O______Description___________________________
C  PIT[N]    Integer     In       Pitches (MIDI # or any numbering system
C                                           that assigns greater numbers
C                                           to faster frequencies.)
C  ACC[N]    Real        Out      Melodic accent
C  N         Integer     In       Length of PIT and ACC
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CONTOUR (PIT, ACC, N)
       IMPLICIT NONE
       INTEGER PIT, N
       REAL ACC
       DIMENSION PIT(N), ACC(N)
       
C           locals
       INTEGER I, CI, CIP1      ! frequency changes at i, i+1
       REAL PI, PIP1            ! approx probabilities of changes at i, i+1
       
C The model assigns values P() to the second and third pitches in a 
C three-tone window.  The accent value is the product of the P() values in 
C two successive windows.  For 2...N-1, values of P() come from values 
C given in the table.  At the beginning and end, single tones are assigned
C P=1 and two tones are assigned [0 1] if the notes are different and
C [.5 .5] if the notes are the same.

       IF (N < 1) RETURN       
       ACC(1) = 1.0                       ! one-tone window
       IF (N < 2) RETURN
       IF (PIT(1) .EQ. PIT(2)) THEN        ! two-tone window
         ACC(1) = 0.5                   ! that is:  ACC(1) = ACC(1) * 0.5
         ACC(2) = 0.5
       ELSE
         ACC(1) = 0.0                   ! that is:  ACC(1) = ACC(1) * 0.0
         ACC(2) = 1.0
       END IF
       IF (N < 3) RETURN
       
       DO I = 2, N-1                         ! three-tone windows
         CI = PIT(I) - PIT(I-1)
         CIP1 = PIT(I+1) - PIT(I)
         IF (CI .EQ. 0 .AND. CIP1 .EQ. 0) THEN
           PI = 0.0
           PIP1 = 0.0
         ELSE IF (CI .NE. 0 .AND. CIP1 .EQ. 0) THEN
           PI = 1.0
           PIP1 = 0.0
         ELSE IF (CI .EQ. 0 .AND. CIP1 .NE. 0) THEN
           PI = 0.0
           PIP1 = 1.0
         ELSE IF (CI > 0 .AND. CIP1 < 0) THEN
           PI = 0.83
           PIP1 = 0.17
         ELSE IF (CI < 0 .AND. CIP1 > 0) THEN
           PI = 0.71
           PIP1 = 0.29
         ELSE IF (CI > 0 .AND. CIP1 > 0) THEN
           PI = 0.33
           PIP1 = 0.67
         ELSE IF (CI < 0 .AND. CIP1 < 0) THEN
           PI = 0.5
           PIP1 = 0.5
         ELSE
           PI = 0.
           PIP1 = 0.
           STOP 'cant happen'
         END IF
         
C Multiply the result for P[i] with the value stored in ACC
         ACC(I) = ACC(I) * PI
         
C Store P[i+1]          
         ACC(I+1) = PIP1
       END DO  ! over three-tone windows
       
C Final two-tone window
       IF (PIT(N-1) .EQ. PIT(N)) THEN
         ACC(N) = ACC(N) * 0.5
       ELSE
         CONTINUE                    ! that is:  ACC(N) = ACC(N) * 1.0
       END IF
       
       RETURN
      END ! of CONTOUR


*11111111111111111111111111111111111111111111111111111111111111111111111      
*                add jitter
*__Name____________Type______In/Out________Description__________________
*  e_t[e]          Real      Both          Event times
*  e_k[e]          Integer   Both          Event kinds
*  e_p1[e]         Integer   Both          Event datum #1
*  e_p2[e]         Integer   Both          Event datum #2 (velocity)
*  e_p3[e]         Real      Both          Event datum #3 (duration)
*  e               Integer   In            Length of event arrays
*11111111111111111111111111111111111111111111111111111111111111111111111
      SUBROUTINE jitter (e_t, e_k, e_p1, e_p2, e_p3, e)
       IMPLICIT NONE
       INTEGER e_k, e_p1, e_p2, e
       REAL e_t, e_p3
       DIMENSION e_t(e), e_k(e), e_p1(e), e_p2(e), e_p3(e)
       
       REAL pitpar, velpar, tim_j, pip
       PARAMETER (pitpar = 0.05,   ! applies to pitch
     &            velpar = 0.10,   ! applies to velocity
     &            tim_j = 0.050,   ! applies to onsets and durations
     &            pip = 0.035)     ! negligible duration  
       INTEGER i
       REAL RANG,                  ! normal random variable
     &      URAND                  ! uniform random variable
       
       DO i = 1, e
         IF (e_k(i) .EQ. 1) THEN  ! jitter times         
           e_t(i) = e_t(i) + tim_j * RANG()
           e_t(i) = MAX(e_t(i), 0.)                     ! don't go negative
         
*                jitter pitches (not too severely)
           IF (URAND() < pitpar) THEN   
             e_p1(i) = e_p1(i) + NINT(RANG())
             e_p1(i) = MIN(MAX(0, e_p1(i)), 127)
           END IF
         
*                velocity and duration
           e_p2(i) = e_p2(i) + NINT(velpar * 16256. * RANG())
           e_p2(i) = MIN(MAX(128, e_p2(i)), 16383)   ! stay in bounds
           e_p3(i) = e_p3(i) + tim_j * RANG()
           e_p3(i) = MAX(pip, e_p3(i))   ! stay above zero
         END IF
           
       END DO  ! next event
       RETURN
      END  ! of jitter


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    duroff - convert durations to NoteOff events
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C__Name_____________Type______In/Out________Description_________________
C  e_t[e]           Real      both          event times
C  e_k[e]           Integer   both          event kinds
C  e_c[e]           Integer   in            event channel
C  e_p1[e]          Integer   in            parameter 1 - pitch
C  e_p2[e]          Integer   in            parameter 2 - velocity
C  e_p3[e]          Real      in            parameter 3 - duration
C  e_v[e]           Integer   in            off velocity
C  e_g[e]           Integer   in            portamento prefix
C  e_f[e]           Integer   in            event source file
C  e_r[e]           Integer   in            event track
C  ind[e]           Integer   neither       permutation vector
C  RSW[3@e]         Real      neither       workspace for sorting
C  JSW[3@e]         Integer   neither       workspace for sorting
C  e                Integer   in            length of event list
C  psi[0:127,0:15]  Integer   work          pointer to sound inception
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE duroff (e_t, e_k, e_c, e_p1, e_p2, e_p3, e_v, e_g, 
     &   e_f, e_r, ind, RSW, JSW, e, psi)
       
       IMPLICIT NONE
       INTEGER e_k, e_c, e_p1, e_p2, e_v, e_g, e_f, e_r, ind, JSW, e, 
     &          psi
       REAL e_t, e_p3, RSW
       DIMENSION e_t(e), e_k(e), e_c(e), e_p1(e), e_p2(e), e_p3(e),
     &  e_v(e), e_g(e), e_f(e), e_r(e), ind(e), RSW(3*e), JSW(3*e),
     &  psi(0:127,0:15)

*                constants
       REAL blip
       PARAMETER (blip = 0.001)

*                local
       INTEGER i, j, 
     &         chan,                     ! channel of an event
     &         indoff, indon,            ! index of last Note-Off, Note-On
     &         kind,                     ! kind of an event
     &         note,                     ! MIDI note # of an event
     &         no                        ! # Note-On events
     
       REAL    told                      ! time of first event
       
       LOGICAL lap                   ! some notes overlap
       
!!!!!!!!!!!!!!!!!!!!!!!!!    Debug
       INTEGER LCOUNT  ! how many times LAP has been true
     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C               begin
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       IF (e < 1) RETURN
       LCOUNT = 0
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      
  10   CONTINUE                ! come here if durations are changed
C  Eliminate zero-duration notes.  
C  Make space by converting bogus events (kind=99) to Note-off.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!       PRINT *, 'in duroff with lcount=', LCOUNT
       DO i = 1, e
         IF (e_k(i) .EQ. 1 .AND. e_p3(i) .LE. 0.) e_k(i) = 0  ! zero-dur
         IF (e_k(i) .EQ. 99) e_k(i) = 0
       END DO

*             sort by kind
       DO i = 1, e
         ind(i) = i
       END DO
       CALL RSORTI (e_k, ind, e, 7, JSW)    ! assuming that n_kind < 128
       CALL RORDER (e_t, ind, e)
       CALL IORDER (e_c, ind, e)
       CALL IORDER (e_p1, ind, e)
       CALL IORDER (e_p2, ind, e)
       CALL RORDER (e_p3, ind, e)
       CALL IORDER (e_v, ind, e)
       CALL IORDER (e_g, ind, e)
       CALL IORDER (e_f, ind, e)
       CALL IORDER (e_r, ind, e)
          
*              count
       IF (e_k(1) < 0) PRINT *, 'duroff: negative event kind!!'

       CALL IBSLE(e_k, e, 0, indoff)     ! index of last Note-Off event
       CALL IBSLE(e_k, e, 1, indon)      ! index of last Note-On event
       no = indon - indoff              ! number of Note-On events       
       IF (no > indoff) PRINT *, 'not enough Note-Off events!'

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             make Note-Off events based on duration of a Note-On
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       DO i = 1, no                  ! index for a note-off
         j = indoff + i               ! index of a note-on
         IF (e_p3(j) .LE. 0.) PRINT *, 'duroff: Zero duration!'
         e_t(i) = e_t(j) + e_p3(j)     ! offtime = ontime + duration
         e_c(i) = e_c(j)               ! same channel
         e_p1(i) = e_p1(j)
         e_p2(i) = e_v(j)               ! off-velocity
         e_p3(i) = 0.             ! null values for duratation,
         e_v(i) = 0               !                 off-velocity,
         e_g(i) = -1              !             and portamento prefix
         e_f(i) = e_f(j)        ! same file
         e_r(i) = e_r(j)        ! same track
       END DO
       
*                mark any leftover Note-Offs as bogus
       DO i = no+1, indoff
         e_k(i) = 99
       END DO       
       
*                   sort by time
       DO i = 1, e 
         ind(i) = i 
       END DO
       CALL MSORTR (e_t, ind, e, RSW, JSW)
       CALL IORDER (e_k, ind, e)
       CALL IORDER (e_c, ind, e)
       CALL IORDER (e_p1, ind, e)
       CALL IORDER (e_p2, ind, e)
       CALL RORDER (e_p3, ind, e)
       CALL IORDER (e_v, ind, e)
       CALL IORDER (e_g, ind, e)
       CALL IORDER (e_f, ind, e)
       CALL IORDER (e_r, ind, e)
       
*               fix negative times
       told = e_t(1)
       IF (told < 0.) THEN
         DO i = 1, e
           e_t(i) = e_t(i) - told
         END DO
        END IF
       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    Lengthen durations up to a repeat NoteOff.
C    Shorten durations up to a repeated NoteOn.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       lap = .FALSE.
       DO j = 0, 15               ! all notes hushed at beginning
         DO i = 0, 127
           psi(i,j) = 0
         END DO
       END DO
       
*        for each event 
       DO i = 1, e
         kind = e_k(i)
         chan = e_c(i)
         note = e_p1(i)
         
         IF (kind .EQ. 0) THEN               
           j = psi(note,chan)
           IF (j > 0) THEN         
             psi(note,chan) = -j                  ! turn off note
           ELSE IF (j < 0) THEN
             j = -j
             e_p3(j) = e_t(i) - e_t(j) - blip          ! extend duration
             IF (e_p3(j) .LE. 0.) e_k(j) = 99
!!!             PRINT *, 'extending duration: ', E_P3(J), '=', E_T(I),
!!!     &         '-', E_T(J), ' for i, j', I, J, ' lap: ', LAP
             lap = .TRUE.
           ELSE
             e_k(i) = 99                         ! bogus note off
             PRINT *, 'duroff:  bogus note off!'
           END IF
           
         ELSE IF (kind .EQ. 1) THEN
           j = psi(note,chan)
           psi(note,chan) = i                 ! turn on note
           IF (j > 0) THEN
             e_p3(j) = e_t(i) - e_t(j) - blip      ! shorten duration
             IF (e_p3(j) .LE. 0.) e_k(j) = 99
!!!             PRINT *, 'shortening duration: ', E_P3(J), '=', E_T(I),
!!!     &         '-', E_T(J), ' for i, j', I, J, ' lap: ', LAP
             lap = .TRUE.
           END IF
         END IF           ! NoteOn or NoteOff ?           
       END DO        ! next event
       
       IF (lap) THEN
         LCOUNT = LCOUNT + 1
!         IF (LCOUNT > 1) PRINT *, 'in DUROFF !!! lcount: ', LCOUNT
         IF (LCOUNT < 999) GO TO 10               ! full do-over
       END IF
       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C            Validate aftertouch
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       DO j = 0, 15               ! all notes hushed at beginning
         DO i = 0, 127
           psi(i,j) = 0
         END DO
       END DO

       DO i = 1, e
         kind = e_k(i)
         chan = e_c(i)
         note = e_p1(i)
         
         IF (kind .EQ. 0) THEN
           psi(note,chan) = 0
         ELSE IF (kind .EQ. 1) THEN
           psi(note,chan) = i
         ELSE IF (kind .EQ. 2) THEN
           IF (psi(note,chan) .EQ. 0) THEN
             e_k(i) = 99    ! bogus aftertouch
             PRINT *, 'duroff:  bogus aftertouch!'
           END IF
         END IF
       END DO

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          Check that all notes are turned off again
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       DO j = 0, 15
         DO i = 0, 127
           IF (psi(i,j) > 0) PRINT *, 'duroff: note (event ', i, ') ',
     &              'left sounding at end of sequence!'
         END DO
       END DO
       RETURN           
      END ! of duroff


*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*       greater than or equal to comparison on arbitrary data 
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      FUNCTION BUFGE (BUF, I, J, M, N)
       IMPLICIT NONE
       INTEGER I, J, M, N
       INTEGER*1 BUF
       DIMENSION BUF(*) 
       LOGICAL BUFGE       
       INTEGER K
       
       BUFGE = .TRUE.
       DO 50 K = 0, MIN(M, N) - 1
         IF (BUF(I+K) - BUF(J+K)) 60, 50, 70
  50   CONTINUE
       IF (M > N) RETURN
       
  60   BUFGE = .FALSE.
  70   RETURN
      END ! of BUFGE


*=======================================================================
*       equality comparison on arbitrary data 
*=======================================================================
      FUNCTION BUFEQ (BUF, I, J, M, N)
       IMPLICIT NONE
       INTEGER I, J, M, N
       INTEGER*1 BUF
       DIMENSION BUF(*) 
       LOGICAL BUFEQ
       INTEGER K
       
       BUFEQ = .FALSE.
       IF (M .NE. N) RETURN
       DO 50 K = 0, N-1
         IF (BUF(I+K) .NE. BUF(J+K)) RETURN
  50   CONTINUE
     
       BUFEQ = .TRUE.
       RETURN
      END  ! of BUFEQ


*-----------------------------------------------------------------------
* The Shell sorting procedure is used to sort the elements of an
* array of arbitrary bytestrings stored in a buffer.  
* The sort is implicit; the answer is returned in permutation vector IP.
* Based on SHELL2 by Alfred Morris, Naval Surface Warfare Center
* Public domain.  This program may be used by any person for any purpose.
*
*_Name_____Type_____In/Out___Description________________________________
* BUF[]    Bytes    in       arbitrary data
* IND[N]   Integer  in       indices to beginings of strings to compare
* LNT[N]   Integer  in       lengths of strings
* IP[N]    Ingeger  out      index array 
* N        Integer  in       number of strings
*-----------------------------------------------------------------------
      SUBROUTINE BUFSRT (BUF, IND, LNT, IP, N)
       IMPLICIT NONE
       INTEGER IND, LNT, IP, N
       INTEGER*1 BUF
       DIMENSION BUF(*), IND(N), LNT(N), IP(N)
       
       INTEGER K(10), IMAX, I, II, KI, J, JMAX, L, LL
       INTEGER INDLL, LNTLL, IPLL             ! temp space
       LOGICAL BUFGE
       DATA K / 1, 4, 13, 40, 121, 
     $  364, 1093, 3280, 9841, 29524 /
     
*             initialize IP
       DO I = 1, N
         IP(I) = I
       END DO     

C             Selection of the increments K(I) = (3**I-1)/2
       IF (N < 2) RETURN
       IMAX = 1
       DO 10 I = 3, 10
         IF (N .LE. K(I)) GO TO 20
         IMAX = IMAX + 1
  10   CONTINUE

C            Stepping through the increments K(IMAX),...,K(1)
  20   I = IMAX
       DO 40 II = 1, IMAX
         KI = K(I)

C             Sorting elements that are KI POSITIONS apart
C               so that A(J).LE.A(J+KI) for J=1,...,N-KI
         JMAX = N - KI
         DO 32 J = 1, JMAX
           L = J
           LL = J + KI
           INDLL = IND(IP(LL))
           LNTLL = LNT(IP(LL))
           IPLL = IP(LL)
  30       IF (BUFGE(BUF,INDLL,IND(IP(L)),LNTLL,LNT(IP(L)))) GO TO 31
           IP(LL) = IP(L)
               LL = L
               L = L - KI
               IF (L > 0) GO TO 30
  31       CONTINUE
           
           IP(LL) = IPLL
  32     CONTINUE

  40   I = I - 1
       RETURN
      END ! of BUFSRT
