%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));