%% Tilting of the dipoles
% Atom Number Density      <= 1250 ppum

%% v_z = 500

% Loop over Lx and Ly values
for Lx = 5:1:10
    for Ly = 5:1:10
        % Initialize OptionsStruct
        OptionsStruct = struct;

        % Assign values to OptionsStruct
        OptionsStruct.NumberOfAtoms           = 101250;     
        OptionsStruct.DipolarPolarAngle       = deg2rad(5);
        OptionsStruct.DipolarAzimuthAngle     = 0;
        OptionsStruct.ScatteringLength        = 75.00;        

        OptionsStruct.TrapFrequencies         = [0, 0, 500];
        OptionsStruct.TrapPotentialType       = 'None';

        OptionsStruct.NumberOfGridPoints      = [128, 128];
        OptionsStruct.Dimensions              = [Lx, Ly];      
        OptionsStruct.TimeStepSize            = 0.005;           % in s
        OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
        OptionsStruct.TimeCutOff              = 1E5;             % in s
        OptionsStruct.EnergyTolerance         = 5E-10;
        OptionsStruct.ResidualTolerance       = 1E-05;
        OptionsStruct.NoiseScaleFactor        = 0.05;

        OptionsStruct.MaxIterations           = 10;      
        OptionsStruct.VariationalWidth        = 1.2;       
        OptionsStruct.WidthLowerBound         = 0.01;
        OptionsStruct.WidthUpperBound         = 12;
        OptionsStruct.WidthCutoff             = 5e-3;

        OptionsStruct.PlotLive                = false;
        OptionsStruct.JobNumber               = (Lx - 5) * (10 - 5 + 1) + (Ly - 4); % Unique JobNumber based on Lx and Ly
        OptionsStruct.RunOnGPU                = true;
        OptionsStruct.SaveData                = true;
        OptionsStruct.SaveDirectory           = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree';
        
        options = Helper.convertstruct2cell(OptionsStruct);
        clear OptionsStruct
    
        solver = VariationalSolver2D.DipolarGas(options{:});
        pot = VariationalSolver2D.Potentials(options{:});
        solver.Potential = pot.trap();
    
        % Run Solver
        [Params, Transf, psi, V, VDk] = solver.run();
    end
end
