From help-octave-request  Sat Feb 20 11:20:12 1993
Subject: octave goodies
From: bardo at ief-paris-sud dot fr (Bardo Muller)
To: help-octave
Cc: bardo at ief-paris-sud dot fr
Date: Sat, 20 Feb 93 18:26:45 +0100

Contents:

1) matlab/octave doc:
2) solver for non-square or ill-conditioned matrices systems + application
3) octave plotter
4) emacs matlab/octave

===============================================================================
1) For those who can't wait for the octave doc:

Take the manual of the old matlab (well, some things may differ,
but most is still ok and you get ideas for naming new routines).
Get it for example from cs.ubc.ca:/mirror1/next/2.0-release/source/matlab.tar.Z
Anyone who got this running on a sun - please e-mail if you succeeded.

===============================================================================
2) For those who can't wait for \ doing non-square matrices:

a\b  means inv(a)*b, or more precise
if a is not square, this is done in a least squares sense.
(explanation by the inventor of matlab)


function x=svdsolve(a,b)
%
% x=svdsolve(a,b)
%
% use it for 
%
%           [A] x] =b] when [A] is a non-square matrix
%                                    or ill-conditioned 
%   
%      source: numerical recipes pp 59
%
if (nargin!=2) 
  error ('svdsolve: need two arguments: matrix a, vector b')
endif
[u,w,v]=svd(a);
for i=1:min(size(w))
  w(i,i)=1./w(i,i);
endfor 
x=v*w'*u'*b;
endfunction


function [a,b]=ratsvd(xdata,ydata,ntop,nbot)
%
% [a,b]=ratsvd(xdata,ydata,ntop,nbot)
%
% Pade approximation of discrete ydata=f(xdata) 
% by a rational function of the form:
%
%      y(x) = sum(1=>ntop)( a(j)*x^(j-1) ) /
%          ( 1 + sum(1=>nbot)( b(j)*x^(j)) )
%
%   xdata,ydata     input data vectors (real or complex)
%   ntop,nbot       number of series terms used in the
%                   numerator and the denominator.
%
%            adapted from ratcof.m originally by
%            H.B.Wilson, Oct 1988
%
ydata=ydata(:);  
xdata=xdata(:); 
m=length(ydata);
nbot=nbot-1;
%
if (nargin==3) 
  nbot=ntop 
endif;
if (nbot+ntop > m)
  disp(sprintf ('\n   ratsvd: system under-determined\n)
  disp(sprintf ('\n   ratsvd: check dmax and consider lowering requested order')
  disp(sprintf ('\n   ratsvd: or supply more data points.\n')
endif
x=ones(m,ntop+nbot);  
x(:,ntop+1)=-ydata.*xdata;
for i=2:ntop 
  x(:,i)=xdata.*x(:,i-1); 
endfor
for i=2:nbot
  x(:,i+ntop)=xdata.*x(:,i+ntop-1); 
endfor
% try ab=inv(x)*ydata, (x needs to be square for octave, so far)
ab=svdsolve(x,ydata); 
a=ab(1:ntop); 
if (nbot == 0)
  b=1;
else
  b=[1,ab(ntop+1:ntop+nbot)']';
endif
endfunction


===============================================================================
3) For those who can't wait for a full implemented plot command:


function pp(plotdata)
%
%  prepare plotfile and take it to a place 
%  where the plotter will find it
% 
errmess='\n   pp: need a single matrix with at least 2 columns\n';
if (nargin != 1)
    error(errmess)
elseif (isstr(plotdata))
    error('\n   pp: can`t plot a string\n')
endif
cols=columns(plotdata);
if (cols<2)
    error(errmess)
endif
save /tmp/oplot plotdata
disp(sprintf('\n.. plot in %d variables ready \n',cols-1))
endfunction

My favourite plot programs:

2-D :	xp, is one program of the MISIM circuit simulation package.
	Get its executable from uwcad.ee.washington.edu
	xp is fast, menu-driven and very complete
	start it by the script below

3-D : 	xprism3 from the khoros distribution
	contours, surface, shading, lightsources..
	its little brother, xprism2 is a decent tool as well.
	both understand directly a file saved by octave (# comment)
	the whole package is quite large, ask your site-administrator
	if it is installed somewhere.

A potential viewer is Xess, a powerful (and commercialized, helas)
X/Motif-based spreadsheet which allows dynamic 2-D and 2-D display
updates by other applications via user-written interfaces.  
Get the non-save version from unix.hensa.ac.uk:/pub/uunet/vendor/ais/


#!/bin/sh
#
# shellfile to start the 2-D plotter xp   
# 
NF=`tail -1 /tmp/oplot | awk "{print NF}"`
case $NF in # clumsy...
2) echo "X Y1" >/tmp/oplot.tmp
   echo "PLOT: 1 Variable found";;
3) echo "X Y1 Y2" >/tmp/oplot.tmp
   echo "PLOT: 2 Variables found";;
4) echo "X Y1 Y2 Y3" >/tmp/oplot.tmp
   echo "PLOT: 3 Variables found";;
5) echo "X Y1 Y2 Y3 Y4" >/tmp/oplot.tmp
   echo "PLOT: 4 Variables found";;
6) echo "X Y1 Y2 Y3 Y4 Y5" >/tmp/oplot.tmp
   echo "PLOT: 5 Variables found";;
7) echo "X Y1 Y2 Y3 Y4 Y5 Y6" >/tmp/oplot.tmp
   echo "PLOT: 6 Variables found";;
*) echo "PLOT: not enough data - exiting..
   exit;;
esac
cat /tmp/oplot | tail +5 >>/tmp/oplot.tmp
/net/speedy/home1/bardo1/xmisim/bin/xp -f /tmp/oplot.tmp 
rm -f /tmp/oplot*

========================================================================
4) For those who now start writing .m files:
Get the matlab-mode for emacs in hub.ucsb.edu:/pub/emacs/matlab-mode.el

========================================================================
What's next ? Who knows of matlab/octave file archives, besides those
of research.att.com:/netlib/matlab and its mirrors ?
Comments to the address below.

Hope to hear from you soon, salut

Bardo
------------
Bardo MULLER                           Phone  : [33] [1] 69 41 78 50    
Institut d'Electronique Fondamentale   Fax    : [33] [1] 60 19 25 93    
Bat. 220 Universite Paris Sud          e-mail : bardo at ief-paris-sud dot fr  
91405 ORSAY CEDEX FRANCE