%%HP:T(3)A(D)F(.);
@ POLYFIT by Stefano Oss
DIR
  MINA
    \<< 1 GETI
      DO ROT ROT GETI 4 ROLL DUP2
        IF >
        THEN SWAP
        END DROP
      UNTIL -64 FS?
      END SWAP DROP
    \>>
  MAXA
    \<< 1 GETI
      DO ROT ROT GETI 4 ROLL DUP2
        IF <
        THEN SWAP
        END DROP
      UNTIL -64 FS?
      END SWAP DROP
    \>>
  ROCO
    \<< \-> A
      \<< A SIZE OBJ\-> DROP \-> n m
        \<< 1 m
          FOR j 1 n
            FOR i A { i j } GET
            NEXT n \->ARRY
          NEXT n EVAL
        \>>
      \>>
    \>>
  PLTF
    \<< ROCO \-> VX VY n
      \<< CLLCD "Plotting Points" 1 DISP
VX MINA SWAP DROP
VX MAXA SWAP DROP
VY MINA SWAP DROP
VY MAXA SWAP DROP \-> xmin xmax ymin ymax
        \<< ERASE
'xmin-10*(xmax-xmin)/131' EVAL
'xmax+10*(xmax-xmin)/131' EVAL
'ymin-10*(ymax-ymin)/65' EVAL
'ymax+10*(ymax-ymin)/65' EVAL \-> xi xa yi ya
          \<< xi xa XRNG yi ya YRNG xi yi R\->C xa ya R\->C BOX
          \>> xmax xmin - 65 / 100 * IP 100 /
              ymax ymin - 32 / 100 * IP 100 / \-> rx ry
          \<< 1 n
            FOR i
              VX i GET
              VY i GET \-> px py
              \<< px rx - py R\->C px rx + py R\->C LINE
                  px py ry - R\->C px py ry + R\->C LINE
              \>>
            NEXT
          \>>
        \>>
      \>>
    \>>
  POLF
    \<< { X EQ pol } PURGE 1 + \-> A N
      \<< A PLTP "Fitting..." 1 DISP A ROCO \-> VX VY n
        \<< 1 N
          FOR i 1 N
            FOR j '\GS(k=1,n,VX(k)^(i+j-2))' EVAL
            NEXT
          NEXT { N N } \->ARRY 1 N
          FOR i '\GS(k=1,n,VY(k)*VX(k)^(i-1))' EVAL
          NEXT N \->ARRY SWAP / \-> C
'pol(X)=\GS(k=1,N,C(k)*X^(k-1))' DEFINE
'pol(X)' EVAL STEQ DRAW
'\v/\GS(k=1,n,(pol(VX(k))-VY(k))^2/(n-N))' EVAL
        \>>
      \>> PICT RCL SWAP \->STR "Chisq. = " SWAP + 1 \->GROB
{ # 5h # 5h } SWAP GOR
2 FIX 'EQ' EVAL \->Q 'Fit' \->TAG SWAP DUP PICT STO \->LCD { } PVIEW
    \>>
END