#! /usr/bin/octave -f
##  -*-octave-*-

## Calculates distance versus release angle for object released from swing
## Copyright (C) 2003  Geraint Paul Bevan

## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



parameters.g = 9.81;		  # gravitational acceleration	[m/s^2]
parameters.h = 0.5;		  # height of pip above ground	[m]
parameters.L = 1.0;		  # distance from pivot to pip	[m]
parameters.theta_max = 2*pi/3;	  # maximum angle of pendulum	[rad]
parameters.theta_rel = pi/2;	  # release angle of pendulum	[rad]
parameters.Vplus = 1.0;		  # extra velocity at release	[m/s]

#############################################################
function [x, y, xdot, ydot] = pendulum (parameters)

  ## position and velocity of a pendulum for a given angle from vertical
  ##
  ## x		[m]	distance along x axis (horizontal)
  ## y		[m]	distance along y axis (vertical)
  ## v_x	[m/s]	velocity along x axis
  ## v_y	[m/s]	velocity along y axis

  g = parameters.g;
  h = parameters.h;
  L = parameters.L;
  theta_max = parameters.theta_max;
  theta_rel = parameters.theta_rel;

  ## tangential velocity of an ideal pendulum
  v_t = sqrt (2*g*L*(cos(theta_rel)-cos(theta_max)));

  ## x and y co-ordinates [m]
  x = L * sin(theta_rel);
  y = L * (1 - cos(theta_rel)) + h;

  ## x and y velocities [m/s]
  xdot = v_t * cos(theta_rel);
  ydot = v_t * sin(theta_rel);
endfunction
#############################################################

function s_x = flight (parameters)
  g = parameters.g;  [xr, yr, xdot, ydot] = pendulum (parameters);
  xdotr = xdot ...
      + parameters.Vplus * parameters.L * cos(parameters.theta_rel);
  ydotr = ydot ...
      + parameters.Vplus * parameters.L * sin(parameters.theta_rel);
  ti = ydotr / parameters.g ...
      + sqrt(ydotr^2 + 2*parameters.g*yr) / parameters.g;
  s_x = xr + xdotr * ti;
endfunction
#############################################################
function plot_results (parameters)
  for j = [1:6]
    i = 0;
    parameters.theta_max = j * pi/12;
    for theta = linspace (-pi/8, 3*pi/4, 100)
      i++;
      parameters.theta_rel = theta;
      data(i, 1  ) = theta;
      if (abs(theta) < parameters.theta_max)
	data(i, 1+j) = flight (parameters);
      else
	data(i, 1+j) = 0;
      endif
      if (data(i, 1+j) <= 0)
	data(i, 1+j) = nan;	# don't plot y=0
      endif
    endfor;
  endfor;
  oneplot;
  gset title "maximum distance versus angle"
  gset xlabel "release angle [degrees]"
  gset ylabel "maximum distance [metres]"
  gset key
  gset grid
  gset xtics 10

  eval(sprintf("gset label \"h = %g m , L = %g m , V+ = %g m/s\" at \
  80,0.25 right ", parameters.h, parameters.L, parameters.Vplus));
  
  gset label "q   =" at 60,3.5 left font "Symbol,20"
  gset label "max" at 61,3.4
  plot (data(:,1)*180/pi, data(:,2), '-;15 deg;',
	data(:,1)*180/pi, data(:,3), '-;30 deg;',
	data(:,1)*180/pi, data(:,4), '-;45 deg;',
	data(:,1)*180/pi, data(:,5), '-;60 deg;',
	data(:,1)*180/pi, data(:,6), '-;75 deg;',
	data(:,1)*180/pi, data(:,7), '-;90 deg;')
  gset output "results.eps";
  gset terminal postscript eps;
  replot;
endfunction
#############################################################
function test_pendulum (parameters)
  n = 101;
  theta = linspace (-parameters.theta_max, parameters.theta_max, n);
  for i = 1 : n
    parameters.theta_rel = theta(i);
    [x(i), y(i), xdot(i), ydot(i)] = pendulum (parameters);
  endfor;

  v_t = sqrt(xdot.^2 .+ ydot.^2);

  data = [x, y, v_t];

  gset autoscale x
  gset autoscale y
  gset nokey
  gset grid xtics
  multiplot(2,2)
  subwindow(1,1)
  gset title "speed versus x and y position"
  gset grid xtics ytics ztics 
  gset xlabel 'x [m]'
  gset ylabel 'y [m]'
  gset zlabel 'speed [m/s]'
  gset view 45 , 45
  gset parametric
  gsplot data title 'pendulum' with points
  gset noparametric
  subwindow(2,1)
  gset title "speed versus angle"
  gset polar
  gset grid polar
  gset size square
  polar_data = [theta', v_t];
  gset xrange [-8:8]
  gset yrange [-8:8]
  gplot polar_data with points
  gset nopolar
  gset grid nopolar
  gset size nosquare
  gset autoscale x
  gset autoscale y
  subwindow(1,2)
  gset title "y position versus x position"
  gset grid xtics ytics 
  gset xlabel 'x position [m]'
  gset ylabel 'y position [m]'
  plot (x, y)
  subwindow(2,2)
  gset title "speed versus x position"
  gset grid xtics ytics 
  gset xlabel 'x position [m]'
  gset ylabel 'speed [m/s]'
  plot (x, v_t);
endfunction
#############################################################
## test_pendulum (parameters); pause;
plot_results (parameters);