clear uIntDecompMatrixEvo uNLDecompMatrix timeGrid thetaMeshGridPlot thetaGridAnimate
clear rGridAnimate rMeshGridPlot pEigenBasis pNLMatrix pNLZeroDecomp pIntMatrixEvo
clear pNLDecompMatrix pNLMatrixSliceTemp vEigenBasis vEigenFncZero vNLDecompMatrix vNLZeroDecomp
clear pEigenFncZero entropyConstXGrid JScaledValues eigenFreqMatrix timePeriod pEigenFncPertTemp
clear vNLMatrix rGrid thetaGrid pConstIntMatrix vBar pIntMatrixIteration
clear smallDivisorsBasicNoZero zeroCorrectionBasic nonResCorrectionBasic
clear pEigenFncIterationLoop1 pEigenFncIterationLoop2
clear pOrigNLMatrix pOrigNLDecompMatrix pOrigNLZeroDecomp uOrigNLDecompMatrix
clear pEigenFncIterationLoop pEigenFncPertTemp2
clear pIntMatrixEvoCorrection



%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %
%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %



%Choose the parameters that will determine the number of eigenfunctions
% that are used. The set of eigenfunctions is a two index set, these
% indicies being "rBasisSize" and "thetaBasisSize"
%The eigenfunctions are indexed as
% phi_(n,m,1) = cos(n*theta)*J_n(lambda_(n,m)*sigmaBar*r)
% phi_(n,m,2) = sin(n*theta)*J_n(lambda_(n,m)*sigmaBar*r)
% where thetaBasis = n for n = 0,1,...,N-1
% and rBasisSize = m  for m = 1,...,M
%Each pair (n,m) corresponds to two different eigenfunctions
% (the cos and sin variant) This is excepting the case of n = 0, m = 0,
% which corresponds to the constant eigenfunction.
% so there will be a total of (N-1)*M + 1 eigenfunctions in the p-basis and
% effectively (N-1)*M eigenfunctions in the u-basis (because we will always
% have that phi_(0,0,1) = c, we know psi_(0,0,1) = grad(phi_(0,0,1)) = 0)
%However, really we lose another M eigenfunctions, as for n=0 the sin
% variant phi_(0,m,2) = 0 for all m, but these still have a row in
% matricies for simplicity



%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %
%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %



%Set whether to run the initialization part of the code that:
% 1. sets all the constants (background pressure, entropy, number of grid
% points, which eigenfunction we are choosing, number of eigenfunctions, etc.)
% 2. Calculates the basis of eigenfunctions for those constants
perform_Initialization = true;


%Set whether to run an iteration method or to just evolve the starting
% guess nonlinearly and return the result
perform_Basic_Iteration = true;


%Set whether or not to improve the initial guess of pBar + alpha*phi_k
% before iteration starts, (this means the O(alpha^2) terms that come
% from (1/2)*alpha^2*D^2f(pBar)[phi_k,phi_k]) and if so whether to still
% evolve the uncorrected initial guess anyway
performInitialGuessCorrection = false;
performUncorrectedEvoAsWell = true;
number_Iterations = 10;
% NOTE: to use the initialization correction function, one must ensure that
%  if we are perturbing the eigenfunction (k1,k2,1) or (k1,k2,2), then
%  2*k1 is within our calculated eigen basis, i.e. that 2*k1 < N


%Set whether the main iteration scripts will prompt the user to choose the
% theta-index and r-index of the fixed k-mode
prompt_User_K_Mode_Indicies = true;


%Set whether the main iteration scripts will prompt the user to choose the
%size of the eigenfunction basis (thetaBasisSize and rBasisSize)
prompt_User_Eigen_Basis_Sizes = true;



%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %
%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if perform_Initialization == true
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%prompt user to choose between an isentropic disk (constant entropy)
% or a disk with a radially piecewise constant entropy profile
%This is done by prompting the user to specify the number of entropy levels
% input of 0 = isentropic case
% input of n = radial profile s(r) with n-1 jumps between n levels

    prompt = ['Specify the number of distinct levels for a radial piecewise constant entropy profile s(r) (1 corresponds to the isentropic case): '];
    radial_Entropy_Levels = input([prompt,newline])


    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    if radial_Entropy_Levels == 1
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

        run("TwoD_Disk_TimeEvo_Isentropic_Main_Iteration_v2.m")

    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    end
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %



    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    if radial_Entropy_Levels > 1
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

        run("TwoD_Disk_TimeEvo_NonIsentropic_Iteration_Main.m")

    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    end
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %
%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %

Error using input
Cannot call INPUT from EVALC.

Error in TwoD_Disk_TimeEvo_Main (line 96)
    radial_Entropy_Levels = input([prompt,newline])
                            ^^^^^^^^^^^^^^^^^^^^^^^

Published with MATLAB® R2025a