Dynamics Notes

Thu Jan 20

     cd ~tpl/matlab
     matlab

     dynamics

Time taken so far 4 hrs.

I'm managing to crash the original program... perhaps no surprise now I've read some of the comments - "QWORK1 mysteriously zeroed by here sometimes", etc.

Fri Jan 21

2 hours today.

Tue Jan 25

• The NAG F04ADF routine "calculates the approximate solution of a set of complex linear equations with multiple right-hand sides, using an LU factorisation with partial pivoting"

  In F04ADF(A, x, B, x,x,x, C...), A and B are the matrices. The answer goes into C. A is overwritten by intermediate values, but dynamics doesn't make use of this, so

     CALL F04ADF( QWORK6, I4T, QWORK5, I4T, I4T, I4T, QWORK7, I4T, WORK2, IFL )
  

  becomes

     QWORK7=QWORK6\QWORK5;
  

  The NAG example provided for the routine solves AX=B where

  A=[ 1, 1+2i, 2+10i; 1+i, 3i, -5+14i; 1+i, 5i, -8+20i] 
  B=[ 1;0;0]
  

  Their solution is

  (10.0000, 1.0000)
  ( 9.0000,-3.0000)
  (-2.0000, 2.0000)
  

  In matlab, x = A\B works giving

  x =
    10.0000 + 1.0000i
     9.0000 - 3.0000i
    -2.0000 + 2.0000i
  

• F01ACF was Withdrawn at Mark 16. Replaced by F01ABF

  Old: CALL F01ACF(N,EPS,A,IA,B,IB,Z,L,IFAIL)
  New: CALL F01ABF(A,IA,N,B,IB,Z,IFAIL)
  

  Given

  A=[ 5 7 6 5; 7 10 8 7; 6 8 10 9; 5 7 9 10];
  

  F01ABF returns

   68
  -41 25
  -17 10  5
   10 -6 -3 2
  

  Matlab's tril(inv(A)) does the same.
• F02AGF was later withdrawn.

  Old: CALL F02AGF(A,IA,N,RR,RI,VR,IVR,VI,IVI,INTGER,IFAIL)
  New: CALL F02EBF('V',N,A,IA,RR,RI,VR,IVR,VI,IVI,WORK,LWORK,IFAIL)
  

  Given

  A =
      0.3500    0.4500   -0.1400   -0.1700
      0.0900    0.0700   -0.5400    0.3500
     -0.4400   -0.3300   -0.0300    0.1700
      0.2500   -0.3200   -0.1300    0.1100
  

  from the NAG example, matlab's

     a=eig(A)
  

  returns the 4 complex eigenvalues. Doing

    [V,D] = eig(A)
  

  returns eigenvectors in V. 3 of the 16 elements have the wrong sign when compared with the NAG document - a typo I presume.

Fri Jan 28

• set m1, c1, c2, k1, k2 all to 10
• set x1 and d/dx1 to 10 in the Initial Conditions
• ran the calculation

• Displacement - Y axis in meters, 0-11. Starts at 10, peaks to 11 at about a sec, then is 0 at about 3.5 secs
• Acceleration - Y axis -40 to 5. Starts at 0. Peaks to 5 at 2secs. After 4 secs flattens to 0
• Velocity - Y axis -10 to 10. Starts at 10. 0 at 1 sec. After 4 secs flattens to 0

I've created a directory for the tcalcs routines and have started converting them.

In that situation I think the original program does

 CALL  TIMCAL(IANAL,INPTS,
     &                NF,NDRFL,
     $                RMASS,
     $                X1DOT,DX,
     o                RTOL,T)

Time today: 3 hours

Mon Jan 31

Tue Feb 1

Wed Feb 9

D02BBF
Withdrawn at Mark 18
Replaced by D02PCF and associated D02P routines

Old: CALL D02BBF(X,XEND,N,Y,TOL,IRELAB,FCN,OUTPUT,W,IFAIL)
New: CALL D02PVF(N,X,Y,XEND,TOL,THRES,2,'usualtask',.FALSE.,
    +            0.0e0,W,14*N,IFAIL)
     ... set XWANT ...
  10 CONTINUE
     CALL D02PCF(FCN,XWANT,X,Y,YP,YMAX,W,IFAIL)
     IF (XWANT.LT.XEND) THEN
       ... reset XWANT ...
       GO TO 10 
     ENDIF

function D02BBF
TSTART=0;
YSTART(1)=0;
YSTART(2)=1;
TEND=2*pi;
NPTS=8;
TINC=(TEND-TSTART)/NPTS;
T=TSTART:TINC:TEND
[t,y]=ode45(@F,T,YSTART)

function YP=F(T,Y)
YP(1)=Y(2);
YP(2)=-Y(1);
YP=YP';

Thu Feb 10

• Displacement - dealt with by grftim.f:DISPLS
• Acceleration- dealt with by grftim.f:ACCNS
• Velocity - dealt with by grftim.f:VELS

      FUNCTION ISTART(NPARAM,NSELC,NPPCOR)
C     NPARAM     The number of the coordinate info required.
C                Numbers refer to coordinates as follows
C                for stage 1
C                                          1 - X1
C                                          2 - X2
C                                          3 - X3
C                                          4 - X1 velocity
C                                          5 - X2 velocity
C                                          6 - X3 velocity
C                                          7 - X1 acceleration
C                                          8 - X2 acceleration
C                                          9 - X3 acceleration
C     NSELC (I)  I=1,9 in stage 1
C                Contains information about whether info about
C                a given coordinate is required.
C                I=1,9 refers to coordinates as in NPARAM
C                NSELC(I) = 1 data reqd
C                NSELC(I) = 0 data not reqd
C    NUMPTS      Number of data points for each coordinate

This program will be harder to convert than previous ones because the fortran maps less cleanly onto matlab.

• There are far more global variables. Many numerical variables have character shadow variables, though from their names this isn't always obvious - e.g. CCIT, CCIOTS and CTOL are the character representations of T (integration time), DELTAT (time step) and RTOL (the tolerance required by the solving routine)
• It uses several (withdrawn) NAG routines. These add clutter to the fortran code.
• Some arrays begin with a 0 index
• There are far more (and more various) function calls. Function-call depths of 5 or so are common, especially with the wrapped NAG routines
• Animation might be a problem

Tue Feb 15

Thu Feb 17

For the first time a non-null-interval calculation runs, giving this output for y in numint.

   10.0000   10.0000
   10.1225    9.5023
   10.2387    9.0141
   10.3489    8.5351
   10.4532    8.0654
   10.8890    5.8521
   11.1927    3.8567
   11.3777    2.0682
   11.4566    0.4754
   11.3683   -1.9354
   11.0379   -3.8139
   10.5226   -5.2265
    9.8711   -6.2336
    9.1236   -6.8939
    8.3182   -7.2619
    7.4852   -7.3892
    6.6494   -7.3211
    5.6304   -7.0220
    4.6663   -6.5516
    3.7768   -5.9697
    2.9751   -5.3232
    2.2093   -4.5890
    1.5566   -3.8620
    1.0125   -3.1698
    0.5715   -2.5318
    0.2637   -2.0275
    0.0204   -1.5778
   -0.1662   -1.1839
   -0.3032   -0.8455
   -0.3833   -0.6091
   -0.4390   -0.4062
   -0.4742   -0.2346
   -0.4920   -0.0918
   -0.4951    0.0372
   -0.4839    0.1368
   -0.4621    0.2109
   -0.4325    0.2630
   -0.4018    0.2934
   -0.3685    0.3114
   -0.3337    0.3191
   -0.2986    0.3184
   -0.2544    0.3078
   -0.2123    0.2892
   -0.1732    0.2653
   -0.1031    0.2062
   -0.0734    0.1742
   -0.0486    0.1436
   -0.0284    0.1151
   -0.0142    0.0926
   -0.0030    0.0724
    0.0056    0.0547
    0.0121    0.0394
    0.0159    0.0286
    0.0186    0.0194
    0.0203    0.0115
    0.0212    0.0050
    0.0214   -0.0009
    0.0210   -0.0054
    0.0201   -0.0088
    0.0189   -0.0111
    0.0176   -0.0125
    0.0162   -0.0134
    0.0147   -0.0138
    0.0132   -0.0138
    0.0113   -0.0134
    0.0095   -0.0127
    0.0078   -0.0117
    0.0062   -0.0105
    0.0047   -0.0091
    0.0033   -0.0077
    0.0022   -0.0064
    0.0013   -0.0051
    0.0007   -0.0041
    0.0002   -0.0032
   -0.0002   -0.0025
   -0.0005   -0.0018
   -0.0007   -0.0013
   -0.0008   -0.0009
   -0.0009   -0.0005
   -0.0009   -0.0002
   -0.0009   -0.0000
   -0.0009    0.0001
   -0.0009    0.0003
   -0.0009    0.0004

Wed Feb 23

More values now put into COORD, though Acceleration has only 2 values associated with it. setfun changed so that more than 1 mass can be used without the program crashing. The following show the results for the 1st mass in the three situations

1 Mass	2 Masses	3 Masses

At least the graphs aren't all the same. There needs to be a way to show the behaviour of the other masses. The original program lets the user selectively display them on the one graph.

I've had a look at the Frequency options. The suboptions are

• Forcing -

  Amplitude	Phase
  X0	X0
  f1	f1

  where the first row is to do with the excitation. The GETFEX fortran routine does this.
• Freq Range - which asks for Start/Final frequencies, and number of points (default 100) The GTFVAL fortran routine does this.
• Calc - the HRMCAL fortran routine does this.
• Graphs - Displacement/Velocity/Acceleration. When the Displacement graph is requested, 2 graphs are shown, both with an x axis running from 0 to 10: The Displacement graph has a 0 to 0.6 y-axis, and the other has a "phase degree" y-axis running from -80 to 0. Both graphs show decaying values with a long tail. I think the displacement graph first rises to a little peak. The FRQGRF fortran routine does this.

Mon Feb 28

• 1 Mass - Displacement and Velocity may be correct. Acceleration has roughly the right range of values.
• 2 Masses - Mass 1 displacement might be right. Mass 2 velocity might be right. The acceleration graphs are wrong.
• 3 Masses - Only the mass 1's displacement graph looks ok.

Thu Mar 3

Wed Mar 9

I'm unsure what the 1st 4 args to timcal do. I think they're something to do with setting one's own x,y points.

• IANAL - integer. Type of analysis performed? 1 Time Domain, 2 Frequency domain, 3 ???, 4 Random/Time, 5 Random/Frequency
• INPTS - array of 30 integers
• NF - integer. Similar to NDOFS?
• NDRFL - array of 30 integers

1 hour.

Thur Mar 10

Well, I've procrastinated all I can. Time to make the single-mass situation produce the correct graphs. Guess I now have to understand how the program's supposed to work so I can fix the Acceleration graph. It would be useful to be able to get some numbers out of the old version rather than just graphs. Alas, the old code fails with

/usr/ccs/bin/ld: Unsatisfied symbols:
   f01acf (code)

2 hours.

Fri Mar 11

• calltimcal CALLS tcalcs/timcal
• tcalcs/timcal CALLS tcalcs/eom1 then tcalcs/vibby
• tcalcs/vibby CALLS tcalcs/numint
• tcalcs/numint CALLS ode45 and tcalcs/out, then displays graphs
• ode45 CALLS FCN
• FCN CALLS setfun
• tcalcs/out CALLS setfun

calltimcal -- timcal -- eom1
                     \_ vibby -- numint -- ode45 -- FCN -- setfun
                                        \_ out   ---------/ 
                                        \_ graphs

Wed Mar 16

I've tried changing masses, etc using the interface (e.g.

• m1=100, m2=1000 and everything else 10
• m1=100, m2=1000, c1=0, c2=100, c3=0 and everything else 10
• m1=100, m2=1000, c1=0, c2=100, c3=0, s2=1000, init displ 1= 100, init velocity 2 =100 and everything else 10

Time for a summary

Thur Mar 17

In hrmcal I don't know what MENOP, MPTS or NPTS are for. Should be ints, so I'll guess MENOP=1; MPTS=0; NPTS=0;

In calc1 I need to determine which of the many calculated values are needed later.

2 hours

Wed Mar 23

If LFP is false then much of the code in calc2 isn't executed. Looking in the fortran code I can't see where LFP is ever set to true

Some matters arising from reading the Instructions for the program "DYNAMICS" documentation

• Features -
  • Points - clicking on the window prints the y coord of the line at the x coord of the mouse click
  • Sessions - can be saved/restored
  • c1, c2.. are "dashpoint rates" and are better labelled as lambda
  • Forcing functions - they can be Fetched (loaded in), Amended (points added, deleted, moved) and Attached (attachment points appear in yellow)
• Changes to the documentation - The instructions can be shortened. Could also note that it's a matlab program so standard methods apply to printing and zooming.

  Printing - use "cuedpr" from the matlab command line to print the current figure in a paper-efficient way. For other hardcopy options use the "Print..." and "Page Setup" menu options. To save images into files use the "Save as" menu option. Zooming - under the "Tools" menu you'll find "Zoom in" and "Zoom out" options. When zoom is on, click the left mouse button in a figure to zoom in on the point under the mouse. Click the right mouse button to zoom out.

Tue Mar 29

My ILAB2 variable in calc2.m is an NxN whereas in the original it's a 1xN (where N in the original might be the number of masses). The contents are correct though - with 1 mass it contains 1 number (namely 1, which in STAGE 1 denotes X1). In calc2.m, WORK1 "contains the negative inverse mass matrix" according to comments. Currently, with a single mass of 10, it's

1 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0

1 hour.

Wed Mar 30

Forcing points can now be removed as well as added.

1 hour.

Fri Apr 1

Here's an updated summary

Thu Apr 21

In calc2.m QWORK1 is (IDIM2*IDIM2*4) which is used as if it were QWORK1(IDIM2*2, IDIM2*2). I'm having trouble with the line

   qcopy2( QWORK1, I4, IDIM2, QWORK7, I4, I4, 1, 1, 1, IDIM3 + 1, IDIM2, I4 );

1 hour.

Tue Apr 26

Tue Apr 27

Tue Apr 28

C     M Hardy uses variable dimensioning of multi- dimensional arrays
C     These must be dimensioned at max and passed in empty through a
C     'buffer layer' where no use is made of them

I've changed qmult4 and qmult5 so that they're closer (to the letter rather than the spirit) to the original fortran code.

2 hours

Tue May 3

I think I'll try to get the fortran program running by removing the graphics. It would help a lot to be able to compare, stage by stage, the fortran and matlab calculations. So in /users0/cstaff/tpl/SRC/dynamics3nographics now are reduced versions of 7 of the fortran files, which I hope will suffice to run the frequency-domain stuff. It compiles with f90 and (unlike the original code) copes with added PRINT statements. I haven't yet worked out how to set all the required values without a user-interface, but I'm progressing.

I've just tried the animation option in the original program. Either it's moving too fast or it's not working. I wonder whether it's still being used? And I'd better try a little animation experiment in matlab to see how jerky/flickery it is. Judge for yourself by running this, which uses (I think) the fastest way to animate

axis([0 10 0 10])
xcoords=[1 2 2 1];
ycoords=[1 1 2 2];
p=patch(xcoords, ycoords,'r')
increment=.02
for i=1:300
  xcoords=xcoords+increment;
  ycoords=ycoords+increment;
  set(p,'XData',xcoords,'YData',ycoords);
  drawnow
end
for i=1:300
  xcoords=xcoords-increment;
  set(p,'XData',xcoords)
  drawnow
end

Wed May 4

• Variables set

  MAXNF=3
  MPTS=500
  RMASS(I)=10
  RSPR(I)=10
  RDMP(I)=10
  NPTS=30
  NDOFS=1
  DX(1)=10
  X1DOT(1)=10
  T=10.0
  DELTAT=0.1
  ICORDS(I)=1
  FFOR(I,J)=.FALSE.
  XARRAY(I,J,K)=0.0
  YARRAY(I,J,K)=0.0
  C     REAL XARRAY(I,J,K) array holding excitation values
  C                I=0,3 0 end wall 1-3 masses
  C                K=3 for excitation
  C                J=1 for amplitude,2 for phase
  XARRAY(0,1,3)=10
  XARRAY(1,1,3)=10
  XARRAY(0,2,3)=10
  XARRAY(1,2,3)=10
  RTOL=0.0001
  F1=0
  F2=100
  NMASS=1
  

• Call tree

main -> setzer
        freqan -> hrmcal -> EOM2 
                            DOHCAL -> CALC2 -> CALC1 -> F01ABF
                                                        NAGCKF
                                                        NAGCFF
                                                        F02AGF
                                               QFORM  
                                               QCOPY1
                                               QIDENT
                                               F04ADF
                                               QCOPY2
                                               QCOPY
                                               QMULT2
                                               QCONVT
                                               QMULT1
                                               QMULT
                                               QMULT4
                                               NAGCKF
                                               NAGCFF
                                               QMULT3
                                               QMULT5
                                               QMULT6

3 hours

Fri Aug 26

Thu Sep 1

The white panels are easier to program than the loose boxes of the "model parameters", but math characters can't be used in them, so they're not used for the model parameters. The screen needs a little tidying up, but I need to concentrate on functionality for a while.

The model parameter callbacks, "Run Time Analysis" and "Save session" buttons now work much as their analogues did before, so do the frequency analysis input boxes. The "Load session" button needs a little more work.

Fri Sep 2

A final decision on whether to use this new version needs to be made around the end of September - the IB Lab runs for the first time on Monday 10 October 2005. We do have the fallback of the existing version ...

Mon Sep 5

Thu Sep 15

Tue Sep 20

• Checking my fortran-to-matlab translation
• Trying to get my non-graphics fortran version to produce the same answers as the installed compiled version
• Suggesting that the matlab code for harmonic analysis be rewritten from scratch

December

January 2006

February 2006

March 2006

                         Dynamics

First some background ... Dynamics is an old fortran program. I started 
converting it to matlab line by line but could never get the freq domain 
stuff working. So I started rewriting the code from first principles. 
This approach will produce a much smaller, faster program that's much 
easier to maintain, but it means that need to understand the code.
I'm aiming to produce a matlab script that given the inputs 
that the old program uses, produces graphs and numerical output. I could
then plug the code into the already-written front-end.

I need help with getting the matlab scripts to produce the output that the 
fortran code does. Ideally of course I'd like the computation code written
for me, but any help on improving/extending would I've done so far would
be useful. Some of the code I've written may confuse more than it helps ...


Now a progress report, with code ... I started first by trying to reproduce 
the graphs that the old program produces.

singleode - this tries to reproduce the 1DOF time/displacement of the
            old program. The shape of the graph is roughly right, but 
            there's less damping than there should be, I think

doubleode - this tries to reproduce some graphs (time and freq) in the 
            2DOF situation. The time ones are more believable than the 
            freq ones.

I then tried to reproduce graph that are in the databook. I had more success
there, but I don't know how to extend these to deal with >1 DOF.

springs  - "Step Response of a linear second-order system
            initially at rest" (p.6 of the data book)

springs2 - "Harmonic response of a linear second-order system -
            Magnitude" (top of p.9 of the data book)

springs3 - "Harmonic response of a linear second-order system -
            Phase" (bottom of p.9 of the data book)

springstest - an attempt to produce the diagram at the top of p.9 of the 
            data book "empirically" (without using the formala on p.8 for 
            the max response
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function singleode
% m1u'' + c1u' + (k1)u = F0sin(wt)

% Using
% u ->u1; u'-> u2

% we get

% u2' = - c1*u2 - ((k1)/m1)*u1  + (F0/m1)*sin(w*t)
% u2(0)=0;
% u2'(0)=1; % a=F/m

m1=10;
c1=10; % damping
k1=10; % stiffness
w=5;   % ???
F0=10; % ???



initial = [10; 10];

[T,Y] = ode23(@yprime2, [0, 10], initial);
figure
plot(T,Y);
%yy=fft(Y)
%size(real(yy))
%figure
%plot(real(yy))

function dy = yprime2(t,y)
w=10; % omega (natural frequency?)

%wn=k/n
wn = sqrt(10/10);
lambda=10;
lambdacrit=2*10*wn;
damping=lambda/lambdacrit;
wd=wn/sqrt(1-2*damping^2);

% This is from WWW notes
% dy= [y(2); -10*y(2)-(wd^2)*y(1)+sin(t)];

% This is from CUED
%  dy= [y(2); -10*y(2)- y(1)+sin(wd*t)];
 dy= [y(2); -y(2)- y(1)+sin(wd*t)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function doubleode


% m1u1'' + c1u1' + (k1+k2)u1 = F0sin(wt)
% m2u2'' + c2u2' - k1u1 = 0


% Using
% u1 ->u1plain; u1'-> u1dash
% u2 ->u2plain; u2'-> u2dash

% we get

% u2' = - c1*u2 - ((k1)/m1)*u1  + (F0/m1)*sin(w*t)
% u2(0)=0;
% u2'(0)=1; % a=F/m

% u2dash' = -c2*u2dash/m2 + u1plain*k1/m2

m1=10;
c1=10; % damping
k1=10; % stiffness
w=5;   % ???
F0=10; % ???



initial = [10; 10; 10; 10];

[T,Y] = ode23(@yprime2, [0, 10], initial);
figure
plot(T,Y);
yy=fft(Y)
size(real(yy))
figure
plot(real(yy))

function dy = yprime2(t,y)
w=10; % omega (natural frequency?)

%wn=k/n
wn = sqrt(10/10);
lambda=10;
lambdacrit=2*10*wn;
damping=lambda/lambdacrit;
wd=wn/sqrt(1-2*damping^2);

dy= [y(2); -10*y(2)-(wd^2)*y(1)+sin(t); y(4) ; -10*y(4)/10 + y(1)*10/10];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function springs
m1=10;
c1=10; % damping
k1=10; % stiffness
f=10; % initial force
x=f/k1;
dampingratio=c1/(2*sqrt(k1*m1));
% Without damping, the natural frequency is w
w=sqrt(k1/m1);

% Y=x/(2*dampingratio*sqrt(1-dampingratio*dampingratio))
natfreq=sqrt(k1/m1)
wdamped=natfreq * sqrt(1-dampingratio*dampingratio)

initial = [0; 0];

for para=[0, 0.1, 0.25, 0.5, 1, 1.5]
     [T,Y] = ode23(@yprime2, [0, 20], initial,[],x,w,para);
     plot(T,Y(:,1)/x);
     hold on
end

function dy = yprime2(t,y,x,w,dampingratio)
dy= [y(2); - 2*dampingratio*y(2)*w - y(1)*w*w + x*w*w];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function springs2
m1=10;
k1=10; % stiffness
natfreq=sqrt(k1/m1);
clf

freqrange=0.01:0.1:2;

for dampingratio=[0.1 0.2 0.3 0.5 1 1.5]
     if dampingratio < 1/sqrt(2) ;
       wmax= natfreq*sqrt(1-2*dampingratio^2);
       responsemax=1/(2*dampingratio*sqrt(1-dampingratio^2));
       disp(sprintf('max response for dampingratio=%.3f',dampingratio))
       disp(sprintf(' is %.3f at freq=%.3f', responsemax, wmax))
       plot(wmax,responsemax,'or');
     end
    element=1;
     for freq=freqrange
          y(element)=realform(natfreq, dampingratio, freq);
          element=element+1;
     end
     plot(freqrange,y)
     hold on
end

function y=realform(omegan, dampingratio, omega) 
omegaratio=omega/omegan;
y= 1/sqrt((1-omegaratio^2)^2 + (2*dampingratio*omegaratio)^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function springs3
m1=10;
k1=10; % stiffness
natfreq=sqrt(k1/m1);
clf

freqrange=0.01:0.1:2;

for dampingratio=[0.1 0.2 0.3 0.5 1 1.5] 
element=1;
for freq=freqrange
     y2(element)=tanphi(natfreq, dampingratio, freq);
     element=element+1;
end
y3=atan(y2)*180/pi;

y3(y3>0) = y3(y3>0)-180;
plot(freqrange/natfreq,y3)
hold on
end

function y=tanphi(omegan, dampingratio, omega) 
     omegaratio=omega/omegan;
     y= (-2*dampingratio*omegaratio)/(1-omegaratio^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function springstest

m1=10;
c1=10; % damping
k1=10; % stiffness
f=10;

dampingratio=c1/(2*sqrt(k1*m1)); % not used
x=f/k1;
w=sqrt(k1/m1);


Y=x/(2*dampingratio*sqrt(1-dampingratio*dampingratio));
natfreq=sqrt(k1/m1);
wdamped=natfreq * sqrt(1-dampingratio*dampingratio);

% sprintf('critical damping rate = %d',2*m1*natfreq)
% sprintf('corresponding dampingratio = %d', 2*m1*natfreq/(2*sqrt(k1*m1)))

trange=0:0.1:18;
freqrange=0.01:0.1:2;
initial = [0; 0];

figure(1)
clf
for dampingratio=[0.1 0.2 0.3 0.5 1 1.5]
  if dampingratio < 1/sqrt(2) ;
     wmax= w*sqrt(1-2*dampingratio^2);
     responsemax=1/(2*dampingratio*sqrt(1-dampingratio^2));

     disp(sprintf('max response for dampingratio=%.3f is %.3f at freq=%.3f',dampingratio, responsemax, wmax))
  end
  
  element=1;
  for freq=freqrange
     [T,Y] = ode23(@yprime2, trange, initial,[],x,w,dampingratio, freq);
     MAXY(element)=max(Y(:,1));
     element=element+1;
  end
  plot(freqrange, MAXY)
  hold on
end


figure(2)
clf
for dampingratio=[0.1 0.2 0.3 0.5 1 1.5]
element=1;
for freq=freqrange
      y(element)=realform(natfreq, dampingratio, freq)
      element=element+1;
end
 plot(freqrange,y)
 hold on
end

  

function dy = yprime2(t,y,x,w,dampingratio, freq)
 dy= [y(2);  - 2*dampingratio*y(2)*w - y(1)*w*w + x*cos(freq*t)];
 
 
 function y=realform(omegan, dampingratio, omega) 
 omegaratio=omega/omegan;
 y= 1/sqrt((1-omegaratio^2)^2 + (2*dampingratio*omegaratio)^2);

July 2006

  cd ~tpl/mechmodelling/canonical
  matlab -r mdof

mass = 0.3*[1,0,0;0,1,0;0,0,1];
damp = 0.1*[2,-1,0;-1,2,-1;0,-1,2];
stiff = 10000*[2,-1,0;-1,2,-1;0,-1,2]; Wmax = 500; % (rad/sec)
N_Sample = 512;
Noise_Factor = 0;
out = 1; %response point
inp = 1; %input point

[H,w,par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor,out,inp)

August 2006

• 4th - It looks as if well be running the old version of dynamics - until Xmas anyway. I've been trying out the new code, comparing the results it gives with those of the HP version. It wants frequencies in terms of "(rad/sec)" rather than "Hz". I assumed that 2*pi rad/sec = 1Hz. I used a simple 1DOF example. In the HP version I set the initial force f1 to 10. I can't see an equivalent in the new version. Here's the code

  mass = 0.3      
  damp = 0.1     
  stiff = 10000 
  N_Sample = 256              
  Noise_Factor = 0;     
  out = 1; %response point       
  inp = 1; %input point    
  Wmax=100*2*pi% (rad/sec = 100Hz???)
  [H,w,par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor, ...
  		       out,inp)
  

  The Hcanonical code is about 70 lines long. One bonus is that the code works out some properties of the system

   182.5742 Natural frequencies
   0.0009   Damping ratio
   3.3333   Amplitude modal constant
   0.0000   Phase modal constant (degrees)
  

  The graphs that it produces might be ok.

   plot(w,180/pi.*angle(H))   
  

  produces a graph with the same shape and y-axis as the HP version's displacement graph, but the x-axis is wrong by about a factor of 5. The other graphs don't closely match the old graphics unless they're inverted etc.

  I then tried a 3-DOF example - all masses 0.3, all dampers 0.1 and all springs 10000. I'm assuming that the "Force input" option selects which mass a force (of unspecified magnitude) will be applied to, and that the "Response" option select which mass to display graphs of. Using guesswork and reverse engineering I've found out how to make the new program produce the same graphs as the old program does (using the imaginary part of the Nyquist graph and the Phase part of the Bode graph, and sometimes taking the abs() values). This only covers the Freq domain stuff (not the Time domain). The "Option FRF" can be set to "Receptance", "Mobility", or "Accelerance" which might just correspond to displacement, velocity and acceleration.

  It's about 80 lines of matlab code - the old dynamics program has over 35000 lines of code.

• 9th - I've been trying to perform the experiment using my code and the new matlab code. Using this code with the new matlab

  mass = 45000;      
  damp = 10000;
  stiff = 4000000; 
  N_Sample = 256;              
  Noise_Factor = 0;     
  out = 1; %response point       
  inp = 1; %input point    
  Wmax=3*2*pi% (rad/sec = 3Hz???)
  [H,w,par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor, ...
  		       out,inp);
  clf
  subplot(4,1,1),plot(w/(2*pi),180/pi.*angle(H))   
  subplot(4,1,2),plot(w/(2*pi),20*real(log10(H)))
  ylabel('mobility');
  M = H.*(j.*w); 
  subplot(4,1,3),plot(w/(2*pi),20*real(log10(M)))
  ylabel('velocity');
  A = H.*-(w.^2);
  subplot(4,1,4),plot(w/(2*pi),20*real(log10(A)))
  ylabel('accelerance');
  

  I get the following results

   1.5005  Natural frequencies
   0.0118   Damping ratio
  

  and the following graphs

  

  which have some similarities with the graphs produced by the old dynamics program. E.g. I think graphs 1 and 2 in the following correspond to graphs 2 and 1 in the above.

  

I think I'll need help to run the transient response part of the lab with the old software.
• 14th - Here's a summary of what the matlab-based versions do
  • Matlab conversion - this is the result of converting from the Fortran/NAG routines. The GUI is complete.
  • Other Matlab code - code written from scratch using CUED's databook, etc. No GUI
  • Matlab public-domain - code from Matlab Central. GUI doesn't let users choose the number of masses (though the code does).

  	Matlab conversion	Other Matlab code	Matlab public-domain
  1 DOF time	Might be working
  2 DOF time	Might be working
  1 DOF freq		Might be working	Might be working
  2 DOF freq			Might be working

  The temptation is to add the "Matlab public-domain" code to the "Matlab conversion" program. That's what ~tpl/bin/newdynamics2 does on the HP teaching system.
• Sept 18th - Some questions for Aylmer

  
  Here's a list of questions regarding the viability of the
  public-domain code that I'd like to use in the 1Bdynamics rewrite.
  
  To use the public-domain routines either
  * copy the code to your machine from my ~tpl/mechmodelling/canonical
    folder on the teaching system.
  or
  * start matlab on the teaching system, then type
      path('~tpl/mechmodelling/canonical',path)
  
  1. Do 
       mdof
  and try the program out. Are the results resonable? Is any of 
  the output the type of output required in the 1Bdynamics experiment?
  
  2. The graphical user interface lacks flexibility. In particular, I
  can't see how to choose 1 single mass. Try the following 1 DOF 
  example using the command line
  
       clf
       mass = 0.3;      
       damp = 0.1 ;     
       stiff = 10000 ; 
       N_Sample = 256 ;              
       Noise_Factor = 0;     
       out = 1; %response point       
       inp = 1; %input point    
       Wmax=100*2*pi; 
       [H,w,par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor, ...
  		       out,inp);
       par
  
  The 4 numbers printed out are
    Natural frequencies
    Damping ratio
    Amplitude modal constant
    Phase modal constant (degrees)
  Are they right?
  
  Can we use the returned data to produce some of the graphs required by the
  1Bdynamics lab? E.g. is
       plot(w,180/pi.*angle(H))   
  meaningful?
  
  3. The following 1DOF example displays more graphs of the results. Are any of
  them meaningful in a 1Bdynamics context?
  
       mass = 45000;      
       damp = 10000;
       stiff = 4000000; 
       N_Sample = 256;              
       Noise_Factor = 0;     
       out = 1; %response point       
       inp = 1; %input point    
       Wmax=3*2*pi%
       [H,w,par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor, ...
  		       out,inp);
       clf
       subplot(4,1,1),plot(w/(2*pi),180/pi.*angle(H))   
       subplot(4,1,2),plot(w/(2*pi),20*real(log10(H)))
       ylabel('mobility');
       M = H.*(j.*w); 
       subplot(4,1,3),plot(w/(2*pi),20*real(log10(M)))
       ylabel('velocity');
       A = H.*-(w.^2);
       subplot(4,1,4),plot(w/(2*pi),20*real(log10(A)))
       ylabel('accelerance');
  

  On 19th Sep I checked whether I'd transplanted the Hcanonical (public domain) code correctly into my matlab code by comparing numerical and graphical output between the original Hcanonical code and my version. All seems ok.
• November
  • Nov 6th - Aylmer and I looked at the programs we have and checked a few results. I'm going to correct the 1DOF Time domain program (but for the initial velocity having the wrong sign, it might be ok). We're hoping for help with the 2DOF Frequency domain code. The program's needed for January.
  • Nov 22nd - Some changes made to the code. I've run 1DOF time-analysis tests, printing out hardcopy from the one and new programs, and all seems ok.
  • Nov 27th - With initial displacement and/or velocities non-zero, all was ok. I hadn't tried the step force before though. Surprizingly (given the amount of fortran and guesses involved), for section 4.2 (1DOF time, with a step force) the new program now produces the same results as the old. For now, simulatingforcing simulates the loading in (and attaching of) a forcing function.
  • Nov 28th - Believable output exists for section 4.1 (1DOF freq). N.B. this new version can currently be run by typing ~tpl/bin/newdynamics on the Teaching System. The start of a replacement for the manual is at http://www-h.eng.cam.ac.uk/help/tpl/programs/dynamics/
  • Nov 29th - 2DOF cases: I've only looked at the options that the lab needs (no velocity or acceleration graphs), and I've not tried all the required values of the damping for the absorber, but all seems ok EXCEPT for the Freq graphs relating to the absorber (used in column 2 of the table on page 6). If the bridge response is correct, I presume the absorber behaviour is being calculated correctly, so maybe I'm displaying the wrong info on the graph.
    For the Time domain, I'm using the conversion from the old program (1000s of lines of matlab code that I don't understand). For the Freq domain I'm using the code I found online (about 50 lines of code that I don't understand)
    Graphics for Time are in tcalcs/numint.m. Graphics for Freq are in callfreqcal2.m.
  • Nov 29th - I've been able run through the handout using the new program. I've checked the new program's numerical output against that produced by the old. They're not exactly the same, but there may be old reasons for that. So I've asked for "dynamics" to be installed. The program could be tidied up. I've also set the defaults so that testing (and doing the lab) is easier. The new program displays more numerical output than the old one did.
• January 2007 - the program was tested on Jan 10th. The following fixes were done on Jan 11-12
  • Removed the "Force" button
  • Removed the message that appeared in the command window when the Point button was clicked on (a pop-up window appears instead)
  • After the help window appears you can now continue using the dynamics program without needing to close the help window
  • The start frequency value was ignored by the code I borrowed to do the freq calculations, so I've changed the code
  • There's more space around the "M2" label
  • Default values for masses, etc are now 0.
  • The number of analysis points chosen by the user was ignored. Now fixed
  • Graphs are bigger
  • The diagram for 2 DOF has been corrected, and y1,y2 dimensions shown
  • The stepforce used in the lab is applied by default.
  • There's a button that displays the step force
  • The frequency window now has a working "Point" button
  • Removed most of the debug info which appeared in the command window
  • Superimposed graphs for the 2DOF cases, and changed the "Points" callbacks
  • Documentation updated
  The program was run with students on Jan 22 after which the following changes were made.
  • After clicking on the Points button the coords are displayed as before, but now the points in question are marked on the graph. The points displayed are those whose x coords are closest to the xcoord of the mouse-click.
  The step function now lasts for 20000 secs rather than 20 (I think). When displaying the step function the graph has '..' rather than '20' as the final label on the x-axes