close

flag=0;
flag1=0;
flag2=0;
flag3=0;
flag4=0;
flag5=0;

while(flag==0)
    var = input('Enter a positive initial production rate (c): ');
    if(var<=0)
        disp('           ')
        disp('ERROR: c must be positive.')
        disp('            ')
    else
        flag=1;
        c=var;
    end    
end
disp('           ')
flag=0;
while(flag==0)
    var = input('Enter a value for T: ');
    if(var<=0)
        disp('      ')
        disp('ERROR: T must be positive.')
        disp('        ')
    else
        flag=1;
        T=var;
    end
end
disp('           ')
disp('One moment please...')
disp('               ')

if(T<=1)
    t1 = linspace(0,T,250);
    u1 = zeros(1,250);
    a1 = T - t1;
    s1 = zeros(1,250) + c;
else
    tswitch = T - 1;
    t3 = linspace(0,tswitch,250);
    u3 = zeros(1,250) + 1;
    a3 = exp(tswitch - t3);
    s3 = c*exp(t3);
    t4 = linspace(tswitch,T,250);
    u4 = zeros(1,250);
    a4 = T - t4;
    s4 = zeros(1,250) + s3(250);
    t1 = [t3 t4(2:250)];
    u1 = [u3 u4(2:250)];
    a1 = [a3 a4(2:250)];
    s1 = [s3 s4(2:250)];
end
while(flag1==0)
    disp('Would you like to plot')
    disp('   1. production rate (x^*)')
    disp('   2. adjoint (lambda)')
    disp('   3. reinvestment percentage (u^*)')
    disp('   4. all three')
    var1=input('Type 1, 2, 3, or 4: ');

    if((var1==1) | (var1==2) | (var1==3) | (var1==4))
        flag1=1;
        disp('         ')
        while(flag2==0)
            disp('Would you like to vary any parameters?')
            disp('1. Yes')
            disp('2. No')
            var2=input('Type 1 or 2: ');
            
            if(var2==1)
                disp('         ')
                flag2=1;
                while(flag3==0)
                    disp('Which parameter would you like to vary?')
                    disp('1. c')
                    disp('2. T')
                    var3=input('Type 1 or 2: ');
                    if(var3==1)
                        disp('        ')
                        while(flag4==0)
                            var4 = input('Enter a second c value: ');
                            if(var4<=0)
                                disp('         ')
                                disp('ERROR: c must be positive')
                                disp('        ')
                            else
                                c2 = var4;
                                flag4=1;
                            end
                        end
                        disp('     ')
                        disp('One moment please...')
                        if(T<=1)
                            t2 = linspace(0,T,250);
                            u2 = zeros(1,250);
                            a2 = T - t2;
                            s2 = zeros(1,250) + c2;
                        else
                            tswitch = T - 1;
                            t3 = linspace(0,tswitch,250);
                            u3 = zeros(1,250) + 1;
                            a3 = exp(tswitch - t3);
                            s3 = c2*exp(t3);
                            t4 = linspace(tswitch,T,250);
                            u4 = zeros(1,250);
                            a4 = T - t4;
                            s4 = zeros(1,250) + s3(250);
                            t2 = [t3 t4(2:250)];
                            u2 = [u3 u4(2:250)];
                            a2 = [a3 a4(2:250)];
                            s2 = [s3 s4(2:250)];
                        end
                        flag3=1;
                    elseif(var3==2)
                        disp('        ')
                        while(flag5==0)
                            var5 = input('Enter a second T value: ');
                            if(var5<=0)
                                disp('         ')
                                disp('ERROR: T must be positive')
                                disp('        ')
                            else
                                T2 = var5;
                                flag5=1;
                            end
                        end
                        disp('     ')
                        disp('One moment please...')
                        if(T2<=1)
                            t2 = linspace(0,T2,250);
                            u2 = zeros(1,250);
                            a2 = T2 - t2;
                            s2 = zeros(1,250) + c;
                        else
                            tswitch = T2 - 1;
                            t3 = linspace(0,tswitch,250);
                            u3 = zeros(1,250) + 1;
                            a3 = exp(tswitch - t3);
                            s3 = c*exp(t3);
                            t4 = linspace(tswitch,T2,250);
                            u4 = zeros(1,250);
                            a4 = T2 - t4;
                            s4 = zeros(1,250) + s3(250);
                            t2 = [t3 t4(2:250)];
                            u2 = [u3 u4(2:250)];
                            a2 = [a3 a4(2:250)];
                            s2 = [s3 s4(2:250)];
                        end
                        flag3=1;
                        
                    else
                        disp('         ')
                        disp('Pardon?')
                        disp('           ')
                    end
                end
            
            elseif(var2==2)
                disp('            ')
                flag2=1;
                if(var1==1)
                    plot(t1,s1)
                    xlabel('Time')
                    ylabel('Production Rate')
                elseif(var1==2)
                    plot(t1,a1)
                    xlabel('Time')
                    ylabel('Adjoint')
                elseif(var1==3)
                    plot(t1,u1)
                    xlabel('Time')
                    ylabel('Reinvestment Percentage')
                    axis([0 T -1 2])
                elseif(var1==4)
                    subplot(3,1,1);plot(t1,s1)
                    subplot(3,1,1);xlabel('Time')
                    subplot(3,1,1);ylabel('Production Rate')
                    subplot(3,1,2);plot(t1,a1)
                    subplot(3,1,2);xlabel('Time')
                    subplot(3,1,2);ylabel('Adjoint')
                    subplot(3,1,3);plot(t1,u1)
                    subplot(3,1,3);xlabel('Time')
                    subplot(3,1,3);ylabel('Reinvestment Percentage')
                    axis([0 T -1 2])
                end
                
            else
                disp('           ')
                disp('Pardon?')
                disp('            ')
            end
        end
        
    else 
        disp('     ')
        disp('Pardon?')
        disp('          ')
    end
end

if(var2==1)
    if(var1==1)
        plot(t1,s1,'b',t2,s2,'g')
        xlabel('Time')
        ylabel('Production Rate')
        legend('First value','Second value',0)
    elseif(var1==2)
        plot(t1,a1,'b',t2,a2,'g')
        xlabel('Time')
        ylabel('Adjoint')
        legend('First value','Second value',0)
    elseif(var1==3)
        plot(t1,u1,'b',t2,u2,'g')
        xlabel('Time')
        ylabel('Reinvestment Percentage')
        if(var3==1)
            axis([0 T -1 2])
        else
            axis([0 max(T,T2) -1 2])
        end
        legend('First value','Second value',0)
    elseif(var1==4)
        subplot(3,1,1);plot(t1,s1,'b',t2,s2,'g')
        subplot(3,1,1);legend('First value','Second value',0)
        subplot(3,1,1);ylabel('Production Rate')
        subplot(3,1,1);xlabel('Time')
        subplot(3,1,2);plot(t1,a1,'b',t2,a2,'g')
        subplot(3,1,2);xlabel('Time')
        subplot(3,1,2);ylabel('Adjoint')
        subplot(3,1,2);legend('First value','Second value',0)
        subplot(3,1,3);plot(t1,u1,'b',t2,u2,'g')
        if(var3==1)
            subplot(3,1,3);axis([0 T -1 2])
        else
            subplot(3,1,3);axis([0 max(T,T2) -1 2])
        end
        subplot(3,1,3);xlabel('Time')
        subplot(3,1,3);ylabel('Reinvestment Percentage')
        
        subplot(3,1,3);legend('First value','Second value',0)
    end
end
