MATLAB function to calculate the energy of double wall carbon nanotubes

function [U dr dt]=DWNT(n1, m1, n2, m2)
% this function uses the chiral vectors of two carbon nanotubes as inputs
%(n1 m1)-inner & (n2 m2)-outer and returns the energy of the tube in kcal
%mol^-1 atom^-1 using an empirical equation found in
%"Double-Wall carbon nanotubes...",
%Bakalis & Zerbetto, Chemical Physics Letters 463 (2008) 139-140
%function written by Tom Moore

%define energy constants
a= 1665.51; %kcal A^6 mol^-1 atom^-1
b= -92.85; %kcal A^3 mol^-1 atom^-1
c= 2.50; %kcal A^1 mol^-1 atom^-1
d= 17.24 ; %kcal A^4 mol^-1 atom^-1
f= -58.56; %kcal A^5 mol^-1 atom^-1

%define graphene lattice constant in A
k= 2.49;

% calculate the chiral angles of both tubes in radians
t1_rad=acos((2*n1+m1)/(2*(n1^2+m1^2+n1*m1)^.5));
t2_rad=acos((2*n2+m2)/(2*(n2^2+m2^2+n2*m2)^.5));

dt=t1_rad-t2_rad %may need to do a <1 check here

% calculate the radius of both tubes
r1=(2.49*(n1^2+m1^2+n1*m1)^.5)/(2*pi);
r2=(2.49*(n2^2+m2^2+n2*m2)^.5)/(2*pi);
dr= r2-r1;

%find the energy of the tube
A= (a/dr^6)+(b/dr^3)+(c/dr);
B=sin(dt)*cos(dt);
C=(d/dr^4)+(f/dr^5);
U=A+B*C