function [M,C,G] = dynamical_properties(q,qd)

q1 = q(1); q2 = q(2); q3 = q(3); 
q1d = qd(1); q2d = qd(2); q3d = qd(3); 

%% phantom premium 1.5 parameters
ma = 0.0202;%new
mc = 0.0249;%new
mbe = 0.2359;%new
mdf = 0.1906;
l1 = 0.215;%new
l2 = 0.170;
l3 = 0.0325;
l5 = -36.8e-3;
l6 = 52.7e-3;
g = 9.81;

%% moment of inertia params
Iaxx = 0.4864e-4;
Iayy = 0.001843e-4;
Iazz = 0.4864e-4;
Icxx = 0.959e-4;
Icyy = 0.959e-4;
Iczz = 0.0051e-4;
Ibexx = 11.09e-4;
Ibeyy = 10.06e-4;
Ibezz = 0.591e-4;
Idfxx = 7.11e-4;
Idfyy = 0.629e-4;
Idfzz = 6.246e-4;
Ibaseyy = 11.87e-4;

%% mass matrix components
m11 = (0.125*(4*Iayy+4*Iazz+8*Ibaseyy+4*Ibeyy+4*Ibezz+4*Icyy+4*Iczz+4*Idfyy+4*Idfzz+4*l1*l1*ma+l2*l2*ma+l1*l1*mc+4*l3*l3*mc)...
    +0.125*(4*Ibeyy-4*Ibezz+4*Icyy-4*Iczz+l1*l1*(4*ma+mc))*cos(2*q2)...
    +0.125*(4*Iayy-4*Iazz+4*Idfyy-4*Idfzz-l2*l2*ma-4*l3*l3*mc)*cos(2*q3)+l1*(l2*ma+l3*mc)*cos(q2)*sin(q3));
m22 = 0.25*(4*(Ibexx+Icxx+l1*l1*ma)+l1*l1*mc);
m23 = -0.5*l1*(l2*ma+l3*mc)*sin(q2-q3);
m32 = m23;
m33 = 0.25*(4*Iaxx+4*Idfxx+l2*l2*ma+4*l3*l3*mc);

M = [m11 0 0;
    0 m22 m23;
    0 m32 m33];


%% coriolis matrix components
c11 = 0.125*(-2*sin(q2)*((4*Ibeyy-4*Ibezz+4*Icyy-4*Iczz+4*l1*l1*ma+l1*l1*mc)*cos(q2)+2*l1*(l2*ma+l3*mc)*sin(q3))*q2d...
    +2*cos(q3)*(2*l1*(l2*ma+l3*mc)*cos(q2)+(-4*Iayy+4*Iazz-4*Idfyy+4*Idfzz+l2*l2*ma+4*l3*l3*mc)*sin(q3))*q3d);
c12 = -0.125*((4*Ibeyy-4*Ibezz+4*Icyy-4*Iczz+l1*l1*(4*ma+mc))*sin(2*q2)+4*l1*(l2*ma+l3*mc)*sin(q2)*sin(q3))*q1d;
c13 = -0.125*(-4*l1*(l2*ma+l3*mc)*cos(q2)*cos(q3)-(-4*Iayy+4*Iazz-4*Idfyy+4*Idfzz+l2*l2*ma+4*l3*l3*mc)*sin(2*q3))*q1d;
c21 = -c12;
c22 = 0;
c23 = 0.5*l1*(l2*ma+l3*mc)*cos(q2-q3)*q3d;
c31 = -c13;
c32 = 0.5*l1*(l2*ma+l3*mc)*cos(q2-q3)*q2d;
c33 = 0;

C = [c11 c12 c13;
    c21 c22 c23;
    c31 c32 c33];

%% gravity vector components
tg1 = 0;
tg2 = g*(l1*(ma+0.5*mc)+l5*mbe)*cos(q2);
tg3 = g*(0.5*l2*ma+l3*mc-l6*mdf)*sin(q3);


G = [0;
    tg2;
    tg3];


