Well,thanks to @Dan.My question seems to be implemented like this:
[n,m]=size(d);%assumes that d is a n x m matrix
[X,Y]=meshgrid(1:n,1:m);%your x-y coordinates
x(:,1)=X(:); % x= first column
x(:,2)=Y(:); % y= second column
f=d(:); % your data f(x,y) (in column vector)
%--- now define the function in terms of x
%--- where you use x(:,1)=X and x(:,2)=Y
fun = @(c,x) c(1)*sin(pi*(x(:,1)+c(2))).*sin(pi*(x(:,2)+c(3))) ./ (25*sin(pi/5*(x(:,1)+c(2))).*sin(pi/5*(x(:,2)+c(3))));
%--- now solve with lsqcurvefit
options=optimset('TolX',1e-6);
c0=[1 0 0];%start-guess here
cc=lsqcurvefit(fun,c0,x,f,[],[],options);
Ifit=fun(cc,x);
Ifit=reshape(Ifit,[n m]);%fitting data reshaped as matrix
surf(X,Y,Ifit);
hold on;
plot3(X, Y, dataArray);