I ended up copying your code and running with some modifications. Looks like it will work a little better if you make the suggestions I made above:
1. Change to TEMP = fft(xt)*smpprd;
2. Remove "2/atdl" from your bottom line (where you determine 'n')
Although, there still does appear to be some dependence on "zeta" on getting a perfect fit (i.e., the bigger zeta gets, the bigger your residuals get). Are you sure the functions you have for the time domain and frequency domain are defined properly?