'---------------------------------------------------------------------
' PROGRAM interpola2                          by G 1 M  1988
'---------------------------------------------------------------------
DIM SHARED xn(50),yn(50)  ' Valori dei nodi
DIM SHARED Nw(50),P(50)   ' Vettori per Newton e Neville
DIM SHARED A(50,50),b(50) ' Tabelle coefficienti, termini noti
DIM SHARED s1(50),s2(50)  ' e costanti per le spline
'---------------------------------------------------------------------
DEF FN f(x)=SIN(x)  
'---------------------------------------------------------------------
' MAIN BODY
'---------------------------------------------------------------------
 CLS  
 fl=0
 LOCATE  4,20 : INPUT "Numero tratti                       ";NumTratti
 LOCATE  6,20 : INPUT "Numero nodi per tratto (Spline=2)   ";NumNodi
 LOCATE  8,20 : INPUT "Inserimento M)anuale F)unzione      ";Modeins$
 
 n=NumNodi-1 : m=NumTratti
 IF Modeins$="m" THEN 
   lin=10 
   FOR j=0 TO m*n
      LOCATE lin,20 : PRINT USING"Nodo ##";j 
      LOCATE lin,40 : INPUT "x ";xn(j)
      LOCATE lin,50 : INPUT "y ";yn(j)
      lin=lin+1 : IF lin>23 THEN CLS : lin=4
   NEXT j   
 END IF  
 IF Modeins$="f" THEN
   LOCATE 10,20 : INPUT "Nodo inferiore     ";NInf
   LOCATE 12,20 : INPUT "Nodo superiore     ";NSup 
   h=(NSup-NInf)/(m*n)
   FOR  j=0 TO m*n
       xn(j)=NInf+h*j
       yn(j)=FN f(xn(j))
   NEXT j
 END IF
 
 CLS
 LOCATE 12,20 : PRINT "a) Newton" 
 LOCATE 14,20 : PRINT "b) Neville"
 LOCATE 16,20 : PRINT "c) Spline naturali" 
 LOCATE 18,20 : PRINT "d) Spline Periodiche"
 LOCATE 20,20 : INPUT "scelta     (a-d)    ";met$ 
 CLS
 LOCATE  4,20 : INPUT "G)rafico  V)alore     ";ModeCalc$
 IF ModeCalc$="g" THEN  
     LOCATE  8,20 : INPUT "Estremo inferiore          ";EInf
     LOCATE 10,20 : INPUT "Estremo superiore          ";ESup
     LOCATE 12,20 : INPUT "Passo grafico              ";PGra
     LOCATE 14,20 : INPUT "Fattore di scala verticale ";fv
     GOSUB traccia
 END IF
 IF ModeCalc$="v" THEN  
      c$=""
      WHILE c$<>"u" 
        LOCATE  8,20 : INPUT "Valore x             ";x
        LOCATE 12,20 : PRINT "Risultato            ";
        GOSUB calcola
        INPUT "'u' per uscire ";c$
      WEND  
 END IF 
    
END 'interpola2      
'---------------------------------------------------------------------
traccia:
  
 CLS  
 IF met$="a" OR met$="b" THEN n=NumNodi-1        
 IF met$="c" THEN fspl=1 : n=NumTratti
 IF met$="d" THEN fspl=0 : n=NumTratti
 fh=600/(ESup-EInf) 
 
 IF met$="a" OR met$="b" THEN                 'newton & neville
   FOR tratto=1 TO NumTratti
     b=(tratto-1)*n
     CALL crocetta(b)
     FOR x=xn(b) TO xn(b+n) STEP PGra
       IF met$="a" THEN
                      CALL newton(x)          'traccia newton
                      CALL punto(x,nwt,1)
                   ELSE                      
                      CALL neville(x)         'traccia neville
                      CALL punto(x,nev,1)
       END IF
       IF Modeins$="f" THEN CALL punto(x,FN f(x),3)
     NEXT x
     fl=0       'ricalcola i coefficienti di newton
   NEXT tratto
   CALL crocetta(NumTratti)
 END IF
 
 IF met$="c" OR met$="d" THEN                 'traccia spline
   CALL disponi
   FOR k=0 TO n-1
     crocetta(k)
     FOR x=xn(k) TO xn(k+1) STEP PGra
       CALL spline(x,k)
       CALL punto(x,spl,1)
       IF Modeins$="f" THEN CALL punto(x,FN f(x),3)
     NEXT x
   NEXT k 
   CALL crocetta(n)                                          
 END IF
 RETURN 'traccia
'--------------------------------------------------------------------- 
 SUB punto(x,y,c%) STATIC
 SHARED EInf,fh,fv
   h=INT(10+(x-EInf)*fh)
   v=INT(100-y*fv)
   PSET (h,v),c%     
 END SUB 'punto  
'---------------------------------------------------------------------
 SUB crocetta(NodoTratto) STATIC
 SHARED xn,yn,EInf,fh,fv
   h=INT(10+(xn(NodoTratto)-EInf)*fh)
   v=INT(100-yn(NodoTratto)*fv)
   LINE (h-2,v)-(h+2,v),2
   LINE (h,v-2)-(h,v+2),2
 END SUB 'crocetta  
'-----------------------------------------------------------------
calcola:
   
 IF met$="a" OR met$="b" THEN n=NumNodi-1        
 IF met$="c" THEN fspl=1 : n=NumTratti
 IF met$="d" THEN fspl=0 : n=NumTratti

 IF met$="a" OR met$="b" THEN         'newton  & neville
   tratto=1 
   WHILE x<xn((tratto-1)*n)
     tratto=tratto+1
   WEND  
   b=(tratto-1)*n
   IF met$="a" THEN
                 CALL newton(x)       'calcola newton
                 PRINT nwt
               ELSE
                 CALL neville(x)      'calcola neville
                 PRINT nev
   END IF
 END IF
 
 IF met$="c" OR met$="d" THEN         'calcola spline
   CALL disponi
   k=0
   WHILE x<xn(k)
     k=k+1
   WEND  
   CALL spline(x,k)
   PRINT spl
 END IF
 RETURN 'calcola
'---------------------------------------------------------------------    
SUB neville(x) STATIC
SHARED xn,yn,n,nev,b
 FOR i=0 TO n
   P(i)=yn(i+b)
 NEXT i
 FOR k=1 TO n
   FOR i=0 TO n-k
     P(i)=( (x-xn(i+k+b))*P(i)-(x-xn(i+b))*P(i+1) )/( xn(i+b)-xn(i+k+b) )
   NEXT i
 NEXT k
 nev=P(0)
END SUB 'neville
'---------------------------------------------------------------------
SUB coeff STATIC
SHARED xn,yn,Nw,n,fl,b
 FOR i=0 TO n 
   Nw(i)=yn(i+b)
 NEXT i
 FOR k=1 TO n
   FOR i=n TO k STEP -1
     Nw(i)=( Nw(i)-Nw(i-1) )/( xn(i+b)-xn(i-k+b) )
   NEXT i
 NEXT k
 fl=1             ' coefficienti calcolati
END SUB 'coeff
'---------------------------------------------------------------------
SUB newton(x) STATIC
SHARED xn,Nw,n,fl,nwt,b
 IF fl=0 THEN CALL coeff
 s=Nw(n)
 FOR k=n-1 TO 0 STEP -1
   s=Nw(k)+(x-xn(k+b))*s
 NEXT k
 nwt=s
END SUB 'newton
'---------------------------------------------------------------------
SUB disponi STATIC
SHARED A,xn,yn,b,s1,s2,n,fspl

  FOR i=0 TO n          ' azzera la tabella dei coefficienti
    FOR j=0 TO n
      A(i,j)=0  
    NEXT j  
  NEXT i                ' calcola la tabella (linee seguenti)
  
  h0=xn(1)-xn(0) : h1=xn(2)-xn(1)
  A(1,1)=2*(h0+h1) : A(1,2)=h1 
  b(1)=6*( (yn(2)-yn(1))/h1-(yn(1)-yn(0))/h0 )
  h0=h1
  FOR i=2 TO n-2
      h1=xn(i+1)-xn(i)
      A(i,i-1)=h0 : A(i,i)=2*(h0+h1) : A(i,i+1)=h1
      b(i)=6*( (yn(i+1)-yn(i))/h1-(yn(i)-yn(i-1))/h0 ) 
      h0=h1
  NEXT i
  h1=xn(n)-xn(n-1)
  A(n-1,n-2)=h0 : A(n-1,n-1)=2*(h0+h1)  
  b(n-1)=6*( (yn(n)-yn(n-1))/h1-(yn(n-1)-yn(n-2))/h0 )
  
  IF fspl=0 THEN    ' aggiunge coefficienti per le spline periodiche
     h0=xn(1)-xn(0) : h1=xn(n)-xn(n-1)          
     A(0,0)=2*(h0+h1) : A(0,1)=h0 : A(0,n-1)=h1
     A(1,0)=h0 : A(n-1,0)=h1
     b(0)=6*( (yn(1)-yn(0))/h0-(yn(n)-yn(n-1))/h1 )
  END IF   
     
  CALL risolvi                 ' risolve il sistema
     
  IF fspl=0 THEN               ' spline periodiche
              s2(n)=s2(0)
            ELSE               ' spline naturali
              s2(0)=0
              s2(n)=0
  END IF
  
  FOR i=0 TO n-1               ' calcola le costanti s1
    hi=xn(i+1)-xn(i)
    s1(i)=( yn(i+1)-yn(i) )/hi-( 2*s2(i)+s2(i+1) )*hi/6
  NEXT i
  
END SUB 'disponi
'---------------------------------------------------------------------  
SUB risolvi STATIC 
SHARED A,b,s2,n                ' risolve il sistema
                              
  FOR k=fspl+1 TO n-1
    FOR i=k TO n-1
      FOR j=k TO n-1
        A(i,j)=A(i,j)-A(i,k-1)*A(k-1,j)/A(k-1,k-1)
      NEXT j
      b(i)=b(i)-A(i,k-1)*b(k-1)/A(k-1,k-1)
    NEXT i
  NEXT k
  s2(n-1)=b(n-1)/A(n-1,n-1)
  FOR i=n-2 TO fspl STEP -1
    s=b(i)
    FOR j=i+1 TO n-1
      s=s-s2(j)*A(i,j)
    NEXT j
    s2(i)=s/A(i,i)
  NEXT i
  
END SUB 'risolvi
'-------------------------------------------------------------------- 
SUB spline(x,k) STATIC
SHARED yn,xn,s2,s1,n,spl

  dx=x-xn(k)
  spl=( s2(k+1)-s2(k) )/6/( xn(k+1)-xn(k) )*dx
  spl=(s2(k)/2+spl)*dx
  spl=(s1(k)+spl)*dx
  spl=yn(k)+spl

END SUB 'spline
'---------------------------------------------------------------------