clear c1 n chi sinterm s_n M I qdot n_max c1_max

dimension = 36;

n=1:dimension; %X axis
c1=1:dimension; %Y axis

chi = zeros(dimension,dimension); %preallocate for speed
sinterm = zeros(dimension,dimension); %preallocate for speed
s_n = zeros(dimension,dimension); %preallocate for speed
qdot = zeros(dimension,dimension); %preallocate for speed

for n_index=1:length(n) %X
    for c1_index=1:length(c1)  %Y
        chi(n_index,c1_index) = pi.*c1(c1_index).*n(n_index)./(c1(c1_index)+n(n_index)./3)+1;  %x=n, y=c1
    end
end

for n_index=1:length(n)  %X
    for c1_index=1:length(c1)  %Y
        sinterm(n_index,c1_index) = sin(chi(n_index,c1_index));
    end
end

for n_index=1:length(n)  %X
    for c1_index=1:length(c1)  %Y
        s_n(n_index,c1_index) = chi(n_index,c1_index);
        for k = 1:length(n_index) %Big Pi 
            s_n(n_index,c1_index)=s_n(n_index,c1_index)*(1-(chi(n_index,c1_index)/(k*pi))^2);
        end
    end
end

for n_index=1:length(n)  %X
    for c1_index=1:length(c1)  %Y
        qdot(n_index,c1_index)=s_n(n_index,c1_index)/sinterm(n_index,c1_index);
    end
end

[M,I] = max(qdot(:));
[n_max, c1_max] = ind2sub(size(qdot),I);  %Maximum values

%waterfall(qdot);