
I need to find the area under a trapezoid for a research-related project.  I was able to find the area under the trapezoid in MATLAB using the code:
function [int] = myquadrature(f,a,b) 
% user-defined quadrature function
% integrate data f from x=a to x=b assuming f is equally spaced over the interval 
% use type 
% determine number of data points
npts = prod(size(f));
nint = npts -1; %number of intervals
if(npts <=1)
  error('need at least two points to integrate')
end;
% set the grid spacing
if(b <=a)
  error('something wrong with the interval, b should be greater than a')
else
  dx = b/real(nint);
end;
npts = prod(size(f));

% trapezoidal rule
% can code in line, hint:  sum of f is sum(f) 
% last value of f is f(end), first value is f(1)
% code below
int=0;
for i=1:(nint)
    %F(i)=dx*((f(i)+f(i+1))/2);
    int=int+dx*((f(i)+f(i+1))/2);
end
%int=sum(F);

Then to call "myquadrature" I did:
% example function call test the user-defined myquadrature function
% setup some data

% velocity profile across a channel
% remember to use ? for help, e.g. ?seq 
x = 0:10:2000; 
% you can access one element of a list of values using brackets
% x(1) is the first x value, x(2), the 2nd, etc.
% if you want the last value, a trick is x(end)  
% the function cos is cosin and mean gives the mean value
% pi is 3.1415, or pi
% another hint, if you want to multiple two series of numbers together
% for example c = a*b where c(1) = a(1)*b(1), c(2) = a(2)*b(2), etc.
% you must tell Matlab you want "element by element" multiplication
% e.g.:    c = a.*b    
% note the "."
%
h = 10.*(cos(((2*pi)/2000)*(x-mean(x)))+1); %bathymetry
u = 1.*(cos(((2*pi)/2000)*(x-mean(x)))+1);  %vertically-averaged cross-transect velocity
plot(x,-h)
% set begin and end points for the integration 
a = x(1);
b = x(end);
% call your quadrature function.  Hint, the answer should be 30000.
f=u.*h;
val =  myquadrature(f,a,b);
fprintf('the solution is %f\n',val);

This is great, I got the expected answer of 30000.
 
NOW THE ISSUE IS, I HAVE NO IDEA HOW THIS CODE TRANSLATES TO R.  Here is what I attempted to do, and with error messages, I can tell i'm doing something wrong:
  myquadrature<-function(f,a,b){
npts=length(f)
nint=npts-1
if(npts<=1)
error('need at least two points to integrate')
end;
if(b<=a)
error('something wrong with the interval, b should be greater than a')
else
dx=b/real(nint)
end;
npts=length(f)
_________(below this line, I cannot code)
int=0
for(i in 1:(npts-1))
sum(f)=((b-a)/(2*length(f)))*(0.5*f[i]+f[i+1]+f[length(f)])}
    %F(i)=dx*((f(i)+f(i+1))/2);
    int=int+dx*((f(i)+f(i+1))/2);
end
%int=sum(F);
 
Thank you and any potential suggestions would be greatly appreciated.
 
Dr. Argese.   		 	   		  
	[[alternative HTML version deleted]]

