WCF Factor

WCF Factor

Computing Tools forEngineering Design 

Problem 1

Problem 2

Problem 3

Write a function that will compute angle between the diagonals of the parallelogram with 3 of 4vertices given by: A(3;-1;2), B(0;4;6), C(5;1;3). The input to the function should be thecoordinates of the vertices and the output is angle in degrees.

Problem 4

Problem 5

Problem 6

Solution 

estimatePgramAngle.m 

function angle = estimatePgramAngle(A, B, C)

% Getting two sides of a parallelogram initially

vec1 = B – A;

vec2 = B – C;

% Calculating angle with the help of dot product

angle = acosd(dot(vec1,vec2)/norm(vec1)/norm(vec2)); 

findVector.m 

function x = findVector(a,b,c)

% To estimate the vector that is perpendicular to a and b we find its cross

% product

x = cross(a,b); 

FunDer.m 

function der = FunDer(fun, x0)

% Derivative of the function at x0 is computed using Two-Point central

% difference formula

h = 0.001;

der = (fun(x0+h)-fun(x0-h))/(2*h);  

HW1_.m 

%% PART 1

% Defining vectors of Temperature T and Wind Speed V to find Wind Chill

% Factor WCF for each combination of T and V

T = -20:5:55;

V = 0:5:55;

% Calculating WCF given T and V

WCF = zeros(length(T),length(V));

for i = 1:length(T)

for j = 1:length(V)

WCF(i,j) = WindChillFactor(T(i),V(j));

end;

end;

% Writing it to the file

fid = fopen(‘output.txt’, ‘w+’);

fprintf(fid,’PART 1\n’);

fprintf(fid,’======\n’);

fprintf(fid,’For T = -20 to 55 and V = 0 to 55 with a separation of 5,\n’);

fprintf(fid,’WCF matrix =\n’);

for i=1:size(WCF, 1)

fprintf(fid, ‘%f ‘, WCF(i,:));

fprintf(fid, ‘\n’);

end

fclose(fid);

%% PART 2

format long % Setting format to long for display in command window

clc         % Clearing command window

% Defining function 1 to estimate its derivative at 0.6

f1 = @(x) x^3*exp(2*x);

x10 = 0.6;

% Derivative of the function at 0.6 is computed using Two-Point central

% difference formula

der1 = FunDer(f1, x10);

% Printing the result to the command window

fprintf(‘Analytic derivative for the first part is 5.02\n’);

fprintf(‘Derivative estimated by the user-defined function is %.3f\n\n’,der1);

% Defining function 1 to estimate its derivative at 0.6

f2 = @(x) 3^x/x^2;

x20 = 2.5;

% Derivative of the function at 2.5 is computed using Two-Point central

% difference formula

der2 = FunDer(f2, x20);

% Printing the result to the command window

fprintf(‘Analytic derivative for the first part is 0.744\n’);

fprintf(‘Derivative estimated by the user-defined function is %.3f\n\n’,der2);

% Writing results to the file

fid = fopen(‘output.txt’, ‘a+’);

fprintf(fid,’\n\nPART 2\n’);

fprintf(fid,’======\n’);

fprintf(fid,’Analytic derivative for the first part is 5.02\n’);

fprintf(fid,’Derivative estimated by the user-defined function is %.3f\n\n’,der1);

fprintf(fid,’Analytic derivative for the first part is 0.744\n’);

fprintf(fid,’Derivative estimated by the user-defined function is %.3f\n\n’,der2);

fclose(fid);

%% PART 3

% Defining vertices A, B and C to calculate the angle between them

A = [3, -1, 2]’;

B = [0, 4, 6]’;

C = [5, 1, 3]’;

% Calculate angle between two pairs of lines AB and BC of parallelogram

% given by the vertices A, B and C.

angle = estimatePgramAngle(A, B, C);

% Printing the result to the command window

fprintf(‘The angle between the diagonals of parallelogram is %.2f\n’,angle);

% Writing results to the file

fid = fopen(‘output.txt’, ‘a+’);

fprintf(fid,’\n\nPART 3\n’);

fprintf(fid,’======\n’);

fprintf(fid,’The angle between the diagonals of parallelogram is %.2f\n’,angle);

fclose(fid);

%% PART 4

% Getting three 3-element vectors input from the user

clear a b c

a(1) = input(‘Enter x coordinate of point a between -4 and 4\n’);

a(2) = input(‘Enter y coordinate of point a between -4 and 4\n’);

a(3) = input(‘Enter z coordinate of point a between -4 and 4\n’);

b(1) = input(‘Enter x coordinate of point b between -4 and 4\n’);

b(2) = input(‘Enter y coordinate of point b between -4 and 4\n’);

b(3) = input(‘Enter z coordinate of point b between -4 and 4\n’);

c(1) = input(‘Enter x coordinate of point c between -4 and 4\n’);

c(2) = input(‘Enter y coordinate of point c between -4 and 4\n’);

c(3) = input(‘Enter z coordinate of point c between -4 and 4\n’);

% Calculating vector x

x = findVector(a,b,c);

% Printing the result to the command window

if sum(isnan(x)) == 0

fprintf(‘The vector x which is perpendicular to a and b is:\n’);

disp(x);

else

fprintf(‘The vector x can not be found\n\n’);

end;

% Writing results to the file

fid = fopen(‘output.txt’, ‘a+’);

fprintf(fid,’\n\nPART 4\n’);

fprintf(fid,’======\n’);

fprintf(fid,’Input vectors were:\n’);

fprintf(fid,’a = [%0.2f %0.2f %0.2f]”\n’,a(1),a(2),a(3));

fprintf(fid,’b = [%0.2f %0.2f %0.2f]”\n’,b(1),b(2),b(3));

fprintf(fid,’c = [%0.2f %0.2f %0.2f]”\n’,c(1),c(2),c(3));

if sum(isnan(x)) == 0

fprintf(fid,’Computed vector x is:\n’);

fprintf(fid,’x = [%0.2f %0.2f]”\n’,x(1),x(2));

else

fprintf(‘The vector x can not be found\n’);

end;

fclose(fid);

%% PART 5

% Loading variables from .dat file

load Coordinates.dat

% Printing the result to the command window

fprintf(‘Function to find if a point lies inside or outside of a triangle\n’);

% Getting point O from the user

pointO(1) = input(‘Enter x coordinate of point O\n’);

pointO(2) = input(‘Enter y coordinate of point O\n’);

% Checking whether the point lies inside or outside of the triangle

ifisInside(pointO,Coordinates)

fprintf(‘The point lies inside the triangle\n\n’);

else

fprintf(‘The point is outside of the triangle\n\n’);

end;

% Writing results to the file

fid = fopen(‘output.txt’, ‘a+’);

fprintf(fid,’\n\nPART 5\n’);

fprintf(fid,’======\n’);

fprintf(fid,’Input point was:\n’);

fprintf(fid,’O = [%0.2f %0.2f]”\n’,pointO(1),pointO(2));

ifisInside(pointO,Coordinates)

fprintf(fid,’The point lies inside the triangle\n\n’);

else

fprintf(fid,’The point is outside of the triangle\n\n’);

end;

fclose(fid);

%% PART 6

% Defining parameters

n = 10;

b = pi/3;

a = 0;

% Defining function

f = @(x) sin(x)+cos(x);

% Calculating integral of the function using trapezoidal rule

resTrap = trap(f,a,b,n);

% Calculating integral of the function using simpson’s rule

resSimp = simp(f,a,b,n);

% Printing the result to the command window

fprintf(‘The integral of the given function estimated using Trapezoidal rule from %.3f to %.3f is %.3f\n’,a,b,resTrap);

fprintf(‘The integral of the given function estimated using Simpson”s rule from %.3f to %.3f is %.3f\n’,a,b,resSimp);

% Writing results to the file

fid = fopen(‘output.txt’, ‘a+’);

fprintf(fid,’\n\nPART 6\n’);

fprintf(fid,’======\n’);

fprintf(fid,’The integral of the given function estimated using Trapezoidal rule from %.3f to %.3f is %.3f\n’,a,b,resTrap);

fprintf(fid,’The integral of the given function estimated using Simpson”s rule from %.3f to %.3f is %.3f\n’,a,b,resSimp);

fclose(fid); 

isInside.m 

function inside = isInside(pointO,vertices)

% Separating three vertices into different variables

A = vertices(1,:);

B = vertices(2,:);

C = vertices(3,:);

% If the point lies inside then the sum of areas of the triangles PAB, PBC,

% PAC is equal to the area of triangle ABC

areaABC = findArea(A, B, C);

areaPAB = findArea(pointO, A, B);

areaPBC = findArea(pointO, B, C);

areaPAC = findArea(pointO, A, C);

ifareaABC == areaPAB + areaPBC + areaPAC

inside = 1;

else

inside = 0;

end;

function area = findArea(A, B, C)

area = abs((A(1)*(B(2)-C(2)) + B(1)*(C(2)-A(2))+ C(1)*(A(2)-B(2)))/2.0); 

simp.m 

function res = simp(f,a,b,n)

% Calculating interval

h = (b-a)/n;

% Initializing vector x

for i = 0:n

x(i+1) = a + i*h;

end;

% Getting the overall sum according to simpson’s formula

pSum = f(x(1))+f(x(end));

for i = 1:2:n-1

pSum = pSum + 4*f(x(i+1));

end;

for i = 2:2:n-2

pSum = pSum + 2*f(x(i+1));

end;

res = pSum*h/3;

 trap.m 

function res = trap(f,a,b,n)

% Calculating interval

h = (b-a)/n;

% Initializing vector x

for i = 0:n

x(i+1) = a + i*h;

end;

% Getting the overall sum according to trapezoidal formula

pSum = f(x(1))+f(x(end));

for i = 1:n-1

pSum = pSum + 2*f(x(i+1));

end;

res = pSum*h/2; 

WindChillFactor.m 

function WCF = WindChillFactor(T,V)

% The function of Wind Chill Factor as defined in Homework 1

WCF = 35.7 + 0.6*T – 35.7*V^0.16 + 0.43*T*V^0.16;