home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / library / libry4b.doc < prev    next >
Text File  |  1989-11-10  |  20KB  |  480 lines
  1. .pa
  2.                  QUICK LIST OF MATHEMATICAL SUBROUTINES
  3.  
  4. BFNLQ.... brute force method for nonlinear simultaneous equations (see BRYDN)
  5. BFNLQD... double precision version of BFNLQ (see BRYDN)
  6. BRYDN.... modified Broyden's method for nonlinear simultaneous equations
  7. BRYDND... double precision version of BRYDN
  8. CONJG.... conjugate gradient method for nonlinear simult. equns. (see BRYDN)
  9. CONJGD... double precision version of CONJG (see BRYDN)
  10. CUBIC.... solve a cubic equation
  11. CUBICD... double precision form of CUBIC
  12. FMIN..... find the minimum of a function within an interval
  13. FMIND.... double precision version of FMIN
  14. FZERO.... find the zero (root) of a function within an interval
  15. FZEROD... double precision version of FZERO
  16. GAUSP.... Gauss-Jordan elimination of simultaneous equations
  17. GAUSPD... double precision form of GAUSP
  18. JLDAY.... convert month/day/year to Julian (annual) day
  19. LLSLC.... linear least-squares with linear constraints
  20. LNREG.... linear regression
  21. MADD..... matrix addition
  22. MADDD.... double precision matrix addition
  23. MCPY..... matrix copy
  24. MCPYD.... double precision matrix copy
  25. MDAY..... convert Julian (annual) day to month/day/year
  26. MINV..... matrix inversion
  27. MINVD.... double precision form of MINV
  28. MPRD..... matrix product (multiplication)
  29. MPRDD.... double precision matrix product (multiplication)
  30. MSUB..... matrix subtraction
  31. MSUBD.... double precision matrix subtraction
  32. MTRA..... matrix transpose
  33. MTRAD.... double precision matrix transpose
  34. MTRS..... in-place transpose of a square matrix
  35. MTRSD.... in-place transpose of a double precision square matrix
  36. NTNLQ.... Newton's method for nonlinear simultaneous equations (see BRYDN)
  37. NTNLQD... double precision version of NTNLQ (see BRYDN)
  38. PROOT.... polynomial root finder
  39. QUAD..... solve a quadratic
  40. QUADD.... double precision form of QUAD
  41. RAND..... random number generator (real)
  42. RK4...... 4th order Runge-Kutta (solution of ordinary differential eqns.)
  43. RK6D..... 6th order Runge-Kutta (solution of ordinary differential eqns.)
  44. SPLNE.... determine coefficients for a cubic spline
  45. SVD...... singular value decomposition using Householder transformations
  46. TDIAG.... tridiagonal matrix solver
  47. TDIAGD... double precision form of TDIAG
  48. TREND.... compute trend of periodic function
  49.  
  50. NOTE: A word of caution about matrix operation subroutines.  All of these
  51.       routines (MADD, MSUB, MCPY, MPRD, MSUB, MTRA, MTRS and their double
  52.       precision  counterparts)  use vector emulations.  Subsequently,  no
  53.       arrays  can  cross  a segment boundary.  See  Chapter  7  for  more
  54.       cautions.
  55. .pa
  56.                            MATHEMATICAL SUBROUTINES
  57.  
  58.  
  59. NAME:     BRYDN
  60. PURPOSE:  modified Broyden's method for nonlinear simultaneous equations
  61. TYPE:     subroutine (far external)
  62. SYNTAX:   CALL BRYDN(USER,XMIN,XMAX,X,F,N,M,WORK,MW,MCALL,IPRT,IER)
  63. IMPUT:    USER (far external) a subroutine that YOU MUST PROVIDE and
  64.           must be called by the sequence CALL USER(X,F)
  65.           XMIN,XMAX (REAL*4) arrays, upper and lower limits on X
  66.           N (INTEGER*2) number of elements in X
  67.           M (INTEGER*2) number of residuals
  68.           W (REAL*4) working space of dimension MW
  69.           MW (INTEGER*2) dimension of working space MW=5N+3M+N*N+N*M
  70.           MCALL (INTEGER*2) maximum function calls (set to zero for unlimited)
  71.           IPRT (INTEGER*2) print select (for IPRT=0 nothing is printed,
  72.           for IPRT=1 only a one line summary at the end is printed, for
  73.           IPRT=2 X, F, and G are printed at each iteration, and for IPRT=3
  74.           the Jacobian is also printed at each iteration)
  75.           need to open file 6 on the PC)
  76. OUTPUT:   X (REAL*4) solution vector
  77.           F (REAL*4) final residual
  78.           IER (INTEGER*2) error indicator
  79.                IER=0  indicates normal return
  80.                IER=1  indicates N<2
  81.                IER=2  indicates M<N
  82.                IER=3  indicates XMIN>XMAX conflict
  83.                IER=4  indicates insufficient work space (increase MW)
  84.           note that the printer output is going to go directly to the
  85.           printer unless you first spool it with by CALLing
  86.           SPOOL('FILE.EXT',IER); and no, the redirect ">" will not work
  87. NOTE:     This is a FAR more difficult problem than most people even
  88.           imagine!  You can compare it to finding the lowest spot on
  89.           the face of the Earth given the rules of the game "Battle
  90.           Ship" where all you can do is call out coordinates and your
  91.           opponent answers with the elevation.  You may find a rut or a
  92.           ditch somewhere; but it's highly unlikely that you will find
  93.           death valley, let alone some trench in the Pacific. Now extend
  94.           this analogy to many-dimensional space... get the point?  So
  95.           don't get too upset if BRYDN can't find the solution you want.
  96.  
  97.           I've tried all sorts of algorithms and this seems to work the
  98.           best for general problems.
  99.  
  100.           If you have 6 equations and 6 unknowns then M=N=6.  If you
  101.           have 12 equations and 6 unknowns then M=12 and N=6.
  102.  
  103.           If all of your equations are linear, by all means don't use BRYDN,
  104.           you want LLSLC.  If you have only one unknown then use FMIN.
  105.  
  106.           For double precision use BRYDND.
  107.  
  108.           Also available are the brute force method (BFNLQ & BFNLQD), the
  109.           conjugate gradient method (CONJG & CONJGD), and Newton's method
  110.           (NTNLQ & NTNLQD). The calling sequence, input, and output are
  111.           all the same as for BRYDN & BRYDND. The only differences is the
  112.           method and the required working space, MW (for BFNLQ, MW=3*N+M;
  113.           for CONJG, MW=8*N+M; and for NTNLQ, MW=6*N+2*M+N*N+N*M).
  114.  
  115.           This is a little tricky, so I have included an example at the end of
  116.           this section.
  117.  
  118.  
  119. NAME:     CUBIC
  120. PURPOSE:  solve a cubic equation
  121. TYPE:     subroutine (far external)
  122. SYNTAX:   CALL SUBROUTINE CUBIC(A1,A2,A3,A4,R1,C1,R2,C2,R3,C3)
  123. INPUT:    A1,A2,A3,A4 (REAL*4) coefficients in cubic equation
  124. OUTPUT:   (R1,C1),(R2,C2),(R3,C3) (REAL*4) roots
  125. NOTE:     the equation has the form A1*X**3+A2*X**2+A3*X+A4=0
  126.           for double precision use CUBICD and don't forget to declare
  127.           A1,A2,A3,A4,R1,C1,R2,C2,R3,C3 all to be REAL*8 and if you
  128.           use constants, don't forget the 1.D0 or whatever.
  129. NOTE:     for double precision use CUBICD
  130.  
  131.  
  132. NAME:     FMIN
  133. PURPOSE:  find the minimum of a function within an interval
  134. TYPE:     REAL*4 function (far external)
  135. SYNTAX:   XMIN=FMIN(F,A,B)
  136. INPUT:    A,B (REAL*4) the search interval
  137.           F (a far external REAL*4 function) that you must supply
  138.           that may be called by FX=F(X)
  139. OUTPUT:   location of minimum
  140. NOTE:     for double precision use FMIND
  141.  
  142.  
  143. NAME:     FZERO
  144. PURPOSE:  find the zero (root) of a function within an interval
  145. TYPE:     REAL*4 function (far external)
  146. SYNTAX:   X0=FZERO(F,A,B)
  147. INPUT:    A,B (REAL*4) interval to look for root in
  148.           F (REAL*4 function) that YOU MUST SUPPLY and must be
  149.           called by the sequence FX=F(X)
  150. OUTPUT:   location of the root
  151. NOTE:     for double precision use FZEROD
  152.  
  153.  
  154. NAME:     GAUSP
  155. PURPOSE:  Gauss-Jordan elimination of simultaneous equations
  156. TYPE:     subroutine (far external)
  157. SYNTAX:   CALL GAUSP(A,B,X,JPIVOT,N,D,C,IER)
  158. INPUT:    A,B (REAL*4) arrays containing the equations to be solved
  159.           having dimensions (N,N) and (N) respectively
  160.           JPIVOT (INTEGER*2) working space of dimension "N"
  161.           N (INTEGER*2) the number of equations
  162. OUTPUT:   X (REAL*4) solution vector of dimension N
  163.           D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
  164.           C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
  165.           measure of the gradedness or condition of the matrix
  166.           IER (INTEGER*2) error indicator (IER=0 is normal)
  167. NOTE:     For double precision use GAUSPD and don't forget to declare
  168.           A,B,X,D,C all to be REAL*8.
  169.           Note that A and B will be destroyed upon return.
  170.           Note that A is in row-order as indicated below
  171.  
  172.                A(1,1)*X(1)+A(1,2)*X(2)+...A(1,N)*X(N)=B(1)
  173.                A(2,1)*X(1)+A(2,2)*X(2)+...A(2,N)*X(N)=B(2)
  174.                A(3,1)*X(1)+A(3,2)*X(2)+...A(3,N)*X(N)=B(3)
  175.                ..........................................
  176.                A(N,1)*X(1)+A(N,2)*X(2)+...A(N,N)*X(N)=B(N)
  177.  
  178.           When you really have to be careful about this is with DATA
  179.           statements as FORTRAN fills arrays in column-order.
  180.  
  181.           Full column pivoting is used.
  182.           Vector emulation is used throughout to give maximum speed.
  183.  
  184. NOTE:     for double precision use GAUSPD
  185.  
  186.  
  187. NAME:     JLDAY
  188. PURPOSE:  convert month/day/year to Julian (annual) day
  189. TYPE:     subroutine (far external)
  190. SYNTAX:   CALL JLDAY(MONTH,IDAY,IYEAR,JDAY)
  191. INPUT:    MONTH,IDAY,IYEAR (INTEGER*2)
  192. OUTPUT:   JDAY (INTEGER*2)
  193. NOTE:     I know that "Julian" day is a misnomer; but that's what
  194.           everyone else calls it.
  195.  
  196.  
  197. NAME:     LLSLC
  198. PURPOSE:  linear least-squares with linear constraints
  199. TYPE:     subroutine (far external)
  200. SYNTAX:   CALL LLSLC(X,Y,A,NE,NC,W,NW,IOPT,IER)
  201. INPUT:    X,Y,A,W (REAL*8) see example below
  202.           NE (INTEGER*2) number of variables (X or A)
  203.           NC (INTEGER*2) number of constraints (may be zero)
  204.           NW (INTEGER*2) working space, NW>=(NE+NC)*(NE+NC+3))
  205.           IOPT (INTEGER*2) option, see example below
  206. OUTPUT:   A (REAL*8)
  207.           IER (INTEGER*2) error indicator
  208.                IER=0   normal
  209.                IER=1   dimension error
  210.                IER=2   singular matrix - no inverse
  211.                IER=3   too few constants+constraints (minimum=2)
  212.                IER=4   too few equations (minimum=1)
  213.                IER=5   too few constraints (minimum=0)
  214.                IER=6   too many constraints (maximum=ne-1)
  215.                IER=7   IOPT not defined
  216.                IER=8   matrices not initialized
  217.                IER=9   too few points for least-squares
  218.                IER=10  constraints added before least-squares points
  219.                IER=11  too many constraints were added
  220.                IER=12  too few constraints were specified
  221.                IER=13  working space too small (increase NW)
  222. NOTE:     This is a little tricky, so I have included an example at the end of
  223.           this section. Sorry, there is no REAL*4 equivalent.
  224.  
  225.  
  226. NAME:     LNREG
  227. PURPOSE:  linear regression
  228. TYPE:     subroutine (far external)
  229. SYNTAX:   CALL LNREG(X,Y,N,A,B,R,IER)
  230. INPUT:    X,Y (REAL*4) data points
  231.           N (INTEGER*2) number of points
  232. OUTPUT:   A,B (REAL*4) slope, intercept
  233.           R (REAL*4) regression coefficient in %
  234.           IER (INTEGER*2) error indicator (IER=0 is normal)
  235. NOTE:     LNREG only fits a straight line of the form Y=A*X+B to a set
  236.           of points.  If you want something more elaborate use LLSLC.
  237.  
  238.  
  239. NAME:     MADD
  240. PURPOSE:  matrix addition
  241. TYPE:     subroutine (far external)
  242. SYNTAX:   CALL MADD(A,B,C,N,M)
  243. INPUT:    A (REAL*4) array of dimension A(N,M)
  244.           B (REAL*4) array of dimension B(N,M)
  245.           N (INTEGER*2) number of rows
  246.           M (INTEGER*2) number of columns
  247. OUTPUT:   C (REAL*4) array of dimension C(N,M)
  248. NOTE:     for double precision use MADDD
  249.  
  250.  
  251. NAME:     MCPY
  252. PURPOSE:  matrix copy
  253. TYPE:     subroutine (far external)
  254. SYNTAX:   CALL MCPY(A,B,N,M)
  255. INPUT:    A (REAL*4) array of dimension A(N,M)
  256.           N (INTEGER*2) number of rows
  257.           M (INTEGER*2) number of columns
  258. OUTPUT:   B (REAL*4) array of dimension B(N,M)
  259. NOTE:     for double precision use MCPYD
  260.  
  261.  
  262. NAME:     MDAY
  263. PURPOSE:  convert Julian (annual) day to month/day/year
  264. TYPE:     subroutine (far external)
  265. SYNTAX:   CALL MDAY(JDAY,IYEAR,MONTH,IDAY)
  266. INPUT:    JDAY,IYEAR (INTEGER*2)
  267. OUTPUT:   MONTH,IDAY (INTEGER*2)
  268. NOTE:     I know that "Julian" day is a misnomer; but that's what
  269.           everyone else calls it.
  270.  
  271.  
  272. NAME:     MINV
  273. PURPOSE:  matrix inversion
  274. TYPE:     subroutine (far external)
  275. SYNTAX:   CALL MINV(A,IPIVOT,JPIVOT,N,D,C,IER)
  276. INPUT:    A (REAL*4) array containing the matrix to be inverted
  277.           IPIVOT,JPIVOT (INTEGER*2) working spaces of dimension "N"
  278.           N (INTEGER*2) rank (the number of equations)
  279. OUTPUT:   A (REAL*4) inverted matrix
  280.           D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
  281.           C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
  282.           measure of the gradedness or condition of the matrix
  283.           IER (INTEGER*2) error indicator (IER=0 is normal)
  284. NOTE:     For double precision use MINVD and don't forget to declare
  285.           A,D,C all to be REAL*8.
  286.           Note that A is inverted "in place."
  287.           Note that A is in row-order (see GAUSP).
  288.  
  289.  
  290. NAME:     MPRD
  291. PURPOSE:  matrix product (multiply)
  292. TYPE:     subroutine (far external)
  293. SYNTAX:   CALL MPRD(A,B,C,N,M,L)
  294. INPUT:    A (REAL*4) array of dimension A(N,M)
  295.           B (REAL*4) array of dimension B(M,L)
  296.           N (INTEGER*2) number of rows in A and C
  297.           M (INTEGER*2) number of columns in A and rows in B
  298.           L (INTEGER*2) number of columns in B and C
  299. OUTPUT:   C (REAL*4) array of dimension C(N,L)
  300. NOTE:     for double precision use MPRDD
  301.  
  302.  
  303. NAME:     MSUB
  304. PURPOSE:  matrix subtraction
  305. TYPE:     subroutine (far external)
  306. SYNTAX:   CALL MSUB(A,B,C,N,M)
  307. INPUT:    A (REAL*4) array of dimension A(N,M)
  308.           B (REAL*4) array of dimension B(N,M)
  309.           N (INTEGER*2) number of rows
  310.           M (INTEGER*2) number of columns
  311. OUTPUT:   C (REAL*4) array of dimension C(N,M)
  312. NOTE:     for double precision use MSUBD
  313.  
  314.  
  315. NAME:     MTRA
  316. PURPOSE:  matrix transpose
  317. TYPE:     subroutine (far external)
  318. SYNTAX:   CALL MTRA(A,B,N,M)
  319. INPUT:    A (REAL*4) array of dimension A(N,M)
  320.           N (INTEGER*2) number of rows in A and columns in B
  321.           M (INTEGER*2) number of columns in A and rows in B
  322. OUTPUT:   B (REAL*4) array of dimension B(M,N)
  323. NOTE:     for double precision use MTRAD
  324.  
  325.  
  326. NAME:     MTRS
  327. PURPOSE:  in-place transpose of a square matrix
  328. TYPE:     subroutine (far external)
  329. SYNTAX:   CALL MTRS(A,N)
  330. INPUT:    A (REAL*4) array of dimension A(N,N)
  331.           N (INTEGER*2) number of rows and columns in A
  332. OUTPUT:   A (REAL*4) array of dimension A(N,N)
  333. NOTE:     for double precision use MTRSD
  334.  
  335.  
  336. NAME:     PROOT
  337. PURPOSE:  polynomial root finder
  338. TYPE:     subroutine (far external)
  339. SYNTAX:   CALL PROOT(P,R,C,N,IER)
  340. INPUT:    P (REAL*8) array containing the polynomial coefficients
  341.           N (INTEGER*2) number of terms in P
  342. OUTPUT:   R (REAL*8) array containing the real part of the roots
  343.           C (REAL*8) array containing the complex part of the roots
  344.           IER (INTEGER*2) error indicator
  345.                IER=0 indicates no error
  346.                IER=1 indicates input error
  347.                IER=2 indicates failure to converge
  348. NOTE:     Bairstow's method is used.
  349.           The polynomial must be of the form
  350.  
  351.                P1+P2*X+P3*X**2+P4*X**3+....=0
  352.  
  353.           Also note that this is a VERY difficult problem for large N,
  354.           so don't be too upset if you get complex roots when you
  355.           thought you should be getting real roots.  That's life in the
  356.           real world.
  357.  
  358.  
  359. NAME:     QUAD
  360. PURPOSE:  solve a quadratic
  361. TYPE:     subroutine (far external)
  362. SYNTAX:   CALL SUBROUTINE QUAD(A1,A2,A3,R1,C1,R2,C2)
  363. INPUT:    A1,A2,A3 (REAL*4) coefficients in quadratic equation
  364. OUTPUT:   (R1,C1),(R2,C2) (REAL*4) roots
  365. NOTE:     the equation has the form A1*X**2+A2*X+A3=0
  366.           for double precision use QUADD and don't forget to declare
  367.           A1,A2,A3,R1,C1,R2,C2 all to be REAL*8 and if you use
  368.           constants, don't forget the 1.D0 or whatever.
  369. NOTE:     for double precision use QUADD
  370.  
  371.  
  372. NAME:     RAND
  373. PURPOSE:  return a random number
  374. TYPE:     subroutine (far external)
  375. SYNTAX:   CALL RAND(X)
  376. OUTPUT:   X (REAL*4) a random number between 0 and 1
  377. NOTE:     this uses the clock for a seed and the "old IBM" algorithm which
  378.           has come under fire as of late, but it works just fine for my
  379.           purposes
  380.  
  381.  
  382. NAME:     RK4
  383. PURPOSE:  4th order Runge-Kutta (solve ordinary differential equations)
  384. TYPE:     subroutine (far external)
  385. SYNTAX:   CALL RK4(USER,X,DX,Y,DY,N,W,V)
  386. INPUT:    X (REAL*4) independent variable
  387.           DX (REAL*4) increment in X (step size)
  388.           Y (REAL*4) array, dependent variables at X
  389.           DY (REAL*4) array, change in Y
  390.           N (INTEGER*2) number of dependent variables (Y)
  391.           W (REAL*4) working space of dimension (N,4)
  392.           V (REAL*4) working space of dimension (N)
  393.           USER (far external) a subroutine that YOU MUST PROVIDE and
  394.           must be called by the sequence CALL USER(X,Y,DY)
  395. OUTPUT:   Y (REAL*4) new values of dependent variable
  396. NOTE:     this isn't too tricky; but I have provided an example that may
  397.           help
  398.  
  399.  
  400. NAME:     RK6D
  401. PURPOSE:  6th order Runge-Kutta (solve ordinary differential equations)
  402. TYPE:     subroutine (far external)
  403. SYNTAX:   CALL RK6D(USER,X,DX,Y,DY,N,W,V)
  404. INPUT:    X (REAL*8) independent variable
  405.           DX (REAL*8) increment in X (step size)
  406.           Y (REAL*8) array, dependent variables at X
  407.           DY (REAL*8) array, change in Y
  408.           N (INTEGER*2) number of dependent variables (Y)
  409.           W (REAL*8) working space of dimension (N,8) <-- note that the
  410.           "8" here is not a misprint - it takes 8 steps to get 6th order
  411.           even though it only takes 4 to get 4-order Runge-Kutta
  412.           V (REAL*4) working space of dimension (N)
  413.           USER (far external) a subroutine that YOU MUST PROVIDE and
  414.           must be called by the sequence CALL USER(X,Y,DY)
  415. OUTPUT:   Y (REAL*8) new values of dependent variable
  416. NOTE:     This is just like RK4 except for the "8" above and that you
  417.           need to declare everything REAL*8 instead of REAL*4.
  418.  
  419.  
  420. NAME:     SPLNE
  421. PURPOSE:  determine coefficients for a cubic spline
  422. TYPE:     subroutine (far external)
  423. SYNTAX:   CALL SPLNE(XI,YI,N,C)
  424. INPUT:    XI,YI (REAL*4) points to be fit
  425.           N (INTEGER*2) number of points
  426. OUTPUT:   C (REAL*4) coefficient array of dimension (3,N)
  427. NOTE:     you need to call this once, then use SEVAL to evaluate spline
  428.  
  429.  
  430. NAME:     SVD
  431. PURPOSE:  singular value decomposition using Householder transformations
  432. TYPE:     subroutine (far external)
  433. SYNTAX:   CALL SVD(A,N,M,S,U,V,W,IER)
  434. INPUT:    A (REAL*8) the matrix (N,M) to be decomposed - row order
  435.           N (INTEGER*2) number of rows in A
  436.           M (INTEGER*2) number of columns in A
  437.           W (REAL*8) working space of dimension MAX0(N,M)
  438. OUTPUT:   S (REAL*8) singular values
  439.           U (REAL*8) left eigenvectors (M,N)
  440.           V (REAL*8) right eigenvectors (N,M)
  441.           IER (INTEGER*2) error indicator (IER=0 is normal)
  442.  
  443.  
  444. NAME:     TDIAG
  445. PURPOSE:  tridiagonal matrix solver
  446. TYPE:     subroutine (far external)
  447. SYNTAX:   CALL TDIAG(D,B,X,W,N,IER)
  448. INPUT:    D (REAL*4) matrix of dimension (N,3)
  449.           B (REAL*4) right-hand-side of dimension (N)
  450.           W (REAL*4) working space of dimension (N,3)
  451.           N (INTEGER*2) number of equations
  452. OUTPUT:   X (REAL*4) solution vector of dimension (N)
  453.           IER (INTEGER*2) error indicator (IER=0 is normal)
  454. NOTE:     the tridiagonal system of equations is defined by  [D][X]=[B]
  455.           [D] should look something like the following
  456.  
  457.                  0  2 -2
  458.                 -1  2 -1
  459.                 -1  2 -1
  460.                  .  .  .
  461.                 -2  2  0
  462.  
  463.           For double precision use TDIAGD and don't forget to declare
  464.           D,B,X,W all to be of type REAL*8.
  465.  
  466.  
  467. NAME:     TREND
  468. PURPOSE:  compute trend of periodic function
  469. TYPE:     subroutine (far external)
  470. SYNTAX:   CALL TREND(X,Y,W,NX,N,C,NC)
  471. INPUT:    X,Y (REAL*4) points to be fit
  472.           W (REAL*4) working space of dimension (NX)
  473.           NX (INTEGER*2) number of points
  474.           N (INTEGER*2) number of known points (N<NX)
  475.           NC (INTEGER*2) number of terms in expansion (try 6)
  476. OUTPUT:   Y (REAL*4) upon return, the N+1,N+2,N+3,...NX values of Y will
  477.           be trended
  478.           C (REAL*4) coefficients in trend series
  479. NOTE:     Chebyshev approximation is used
  480.