as i see, you are taking the positive values of m. In other words you ignor the -M values in line code number 41 {for ii = 1 : col
mm = m(ii);}
because u let m=-M:M; your calculations are for the positive values of m (from 1 to M only in line code 41 to 49).
-----------------------------------------
for ii = 1 : col
mm = m(ii);
yy = cos_theta_m(mm, lambda, length, samples);
zz=z_m(mm,d,samples);
sum = sum + (bm(jj)*exp(-j*k*zz*yy));
jj = jj + 1;
AF(kk) = sum;
kk = kk + 1;
end
------------------------------------------
correct me if i am in wrong idea