%Matthew Miller
%1/16/22
%Calculates max qdot for a range of n and c1 values
%to do more I'd have to optimize a time-dependent progression
%I may do it later

%pdot = q1*q2*q
%qdot = m3Factor * s_n(chi)/sin(chi)
%s_n(chi) = chi * PI(k=1:n)(1-(chi/(k*pi))^2)

%==========================================================================
%start fresh

clear c1 n chi sinterm s_n M I qdot n_max c1_max

%==========================================================================
%human entry

dimension = 36;

%not sure what to do with these yet
q1=1;
q2=1;
c1=1;
c2=1;
c3=1;

m1=0;  %Milestone upgrade 1 level (0, 1, 2, 3, or 4)
m2=0;  %Milestone upgrade 2 level (0 or 1)
m3=1;  %Milestone upgrade 3 level (0, 1, 2, or 3)

%==========================================================================
%prelim calculations

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

m1Factor = q1^(1+0.01*m1); %precalculate
m2Factor = c3^m2; %precalculate
m3Factor = 3^m3; %precalculate

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

%==========================================================================
%generate variables
%Normally I'd do all these with functions, but that's harder to read

%chi
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)./m3Factor)+1;  %x=n, y=c1
    end
end

%sin(chi)
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

%s_n(chi)
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

%qdot(chi)
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

%==========================================================================
%find max qdot

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

%waterfall(qdot(1:10,1:10));