; ------------- GAUSS Integrals --------------
;               with five points
; -------- guillaume.tello@orange.fr ---------


; procedures

#label A set_low
#label B set_high
#label C integral
#label E function

; flags

#flag 3 second_part

; variables

#reg 0 incdec_I
#reg 1 J
#reg 2 x
#reg 3 Integ
#reg 4 n
#reg 5 low
#reg 6 high-low
#reg 7 h
#reg 8 I

#array coeff 9 14

;--------------------------------------
; initialized data saved in GAUSS_5.DAT
;--------------------------------------

#register GAUSS_5.DAT

#data

; this table is used from bot to top (3 points)
;      and then from top to bottom (2 more points)

coeff { 0.04691007703067    ; dx1               dx6
        0.118463442528      ; weight of f(x1)   weight of f(x5)
        0.1838552679164     ; dx2               dx5
        0.2393143352496     ; weight of f(x2)   weight of f(x4)
        0.2692346550528     ; dx3               dx4
        0.28444444444444    ; weight of f(x3)
      }

#text

;----------------------------------------
; a A to set the lower bound for integral
;----------------------------------------

LBL set_low
    STO low
RTN

;-----------------------------------------
; b B to set the higher bound for integral
;-----------------------------------------

LBL set_high
    - RCL low =
        STO high-low        ; b-a
RTN

;----------------------------------------------
; n C to calculate integral from a to b of F(x)
;        with n intervals
;----------------------------------------------

LBL integral
    STO n
    
    h={ 1/X * @high-low = }             ; h = (b-a)/n
    x={ @low }                          ; x=a to start
    Integ={ 0 }                         ; integral to zero
    
    FOR n
    
        I={ 8 }
        incdec_I={ $28 }                ; 28 to increment I
        INV STF second_part             ; not second part!
        J={ 3 }                         ; 3 points to calculate in first part
        
        looping:
        
        FOR J
            OP* incdec_I
            x+{ @coeff(I) * @h = }      ; next X
            OP* incdec_I
            Integ+{ @x SBR f(x) * @h * @coeff(I) = }    ; add weighted f(x) to integral
        NEXT
        
        IF(off) second_part             ; if first part
            incdec_I={ $38 }            ; then 38 to decrement I
            J={ 2 }                     ; two points in the second part
            STF second_part             ; put the flag
            GTO looping                 ; and start the loop again
        ENDIF
        
        x+{ RCL 9 * @h = }              ; if second part done, go to next interval
        
    NEXT
    
    RCL Integ                           ; display integral
    
RTN

;------------------------------
; function of x
; X at the display or in Reg 02
;------------------------------

LBL function
f(x):
    EXP                                 ; example! X at the display or in Reg 02.
RTN

#end