Weierstrass Sine Product, qdot simulator.

%Matthew Miller

%v_4

%2/8/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 all

%==========================================================================

%human entry

%set number of purchases

steps = 1000;

%since the maximum number of purchases of each variable is all of the

%purchases, these can be set to steps

nDimension = steps;

c1Dimension = steps;

%not sure what to do with these yet

q1=1;

q2=1;

c2=1;

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

m2=1; %Milestone upgrade 2 level (0 or 1)

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

%==========================================================================

%prelim calculations

nIndex=1:nDimension; %X axis

c1Index=1:c1Dimension; %Y axis

m1Factor = q1^(1+0.01*m1); %precalculate

m2Factor = c2^m2; %precalculate

m3Factor = 3^m3; %precalculate

chi = zeros(nDimension,c1Dimension); %preallocate for speed

sinterm = zeros(nDimension,c1Dimension); %preallocate for speed

s_n = zeros(nDimension,c1Dimension); %preallocate for speed

qdot = zeros(nDimension,c1Dimension); %preallocate for speed

%==========================================================================

%generate variables

%Normally I’d do all these with functions, but that’s harder to read

%blah blah blah arrays begin at 1 so I’m ignoring c1(0)=0

stepLength = 50;

basePower = 1;

offset = 1;

power=basePower;

c1steps=1;

c1(1)=offset;

for index = 2:c1Dimension

c1(index) = c1(index-1)+power;

c1steps=c1steps+1;

if c1steps > stepLength

power=power*2;

c1steps=1;

end

end

%chi

for n_index=1:nDimension %X

for c1_index=1:c1Dimension %Y

chi(n_index,c1_index) = pi.*c1(c1_index).*nIndex(n_index)./(c1(c1_index)+nIndex(n_index)/m3Factor)+1; %x=n, y=c1

end

end

%sin(chi)

for n_index=1:nDimension %X

for c1_index=1:c1Dimension %Y

sinterm(n_index,c1_index) = sin(chi(n_index,c1_index));

end

end

%s_n(chi)

for n_index=1:nDimension %X

for c1_index=1:c1Dimension %Y

s_n(n_index,c1_index) = chi(n_index,c1_index);

s_nk(1)=chi(n_index,c1_index)*(1-(chi(n_index,c1_index)/(1*pi))^2);

if n_index>1

for k = 2:n_index %Big Pi

s_nk(k)=s_nk(k-1)*(1-(chi(n_index,c1_index)/(k*pi))^2);

end

end

s_n(n_index,c1_index)=s_nk(end);

end

end

%qdot(chi)

for n_index=1:nDimension %X

for c1_index=1:c1Dimension %Y

qdot(n_index,c1_index)=m2Factor*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

qdot(n_max,c1_max);

%========================================================================================

%plotter

%========================================================================================

clear rows cols

maxSteps = steps;

minSteps = 1;

for stepsIndex=minSteps:maxSteps

nRange=stepsIndex;

c1Range=stepsIndex;

C=max(qdot(1:nRange,1:c1Range));

D=max(C);

[rows(stepsIndex), cols(stepsIndex)] = find(qdot(1:nRange, 1:c1Range) == D);

end

figure(1)

plot(rows, cols)

title(‘Peak Qdot’)

xlabel(‘n’)

ylabel(‘c1’)