%% Tilting of the dipoles

ppum                                      = 1250; % Atom Number Density in per micrometers

%% theta = 0

LatticeSpacing                            = 1.0:0.05:4.0;
totalIterations                           = numel(LatticeSpacing);

% create a local cluster object
cluster = parcluster('local');

% get the number of dedicated cores from environment
nprocs = str2num(getenv('SLURM_NPROCS'));

% you may explicitly set the JobStorageLocation to the tmp directory that is unique to each cluster job (and is on local, fast scratch)
parpool_tmpdir = [getenv('TMP'),'/.matlab/local_cluster_jobs/slurm_jobID_',getenv('SLURM_JOB_ID')];
mkdir(parpool_tmpdir);
cluster.JobStorageLocation = parpool_tmpdir;

% start the parallel pool
parpool(cluster, nprocs)

% Parallel loop over all combinations of i, j
parfor (k = 1:totalIterations, cluster)
    
    a                                     = LatticeSpacing(k);
    lx                                    = a;
    ly                                    = sqrt(3)*a;
    
    % Initialize OptionsStruct
    OptionsStruct                         = struct;

    % Assign values to OptionsStruct
    OptionsStruct.NumberOfAtoms           = ppum * (lx*ly);     
    OptionsStruct.DipolarPolarAngle       = deg2rad(0);
    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            = 5E-4;            % in s
    OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
    OptionsStruct.TimeCutOff              = 2E6;             % in s
    OptionsStruct.EnergyTolerance         = 5E-10;
    OptionsStruct.ResidualTolerance       = 1E-05;
    OptionsStruct.NoiseScaleFactor        = 0.05;
    OptionsStruct.BiasWithAnsatz          = true;
    OptionsStruct.Ansatz                  = 'triangular';
    OptionsStruct.IncludeDDICutOff        = false;
    
    OptionsStruct.MaxIterations           = 10;      
    OptionsStruct.VariationalWidth        = 1.10;       
    OptionsStruct.WidthLowerBound         = 0.01;
    OptionsStruct.WidthUpperBound         = 12;
    OptionsStruct.WidthCutoff             = 1e-2;

    OptionsStruct.PlotLive                = false;
    OptionsStruct.JobNumber               = k;
    OptionsStruct.RunOnGPU                = true;
    OptionsStruct.SaveData                = true;
    OptionsStruct.SaveDirectory           = sprintf('./Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree%i', round(rad2deg(OptionsStruct.DipolarPolarAngle)));

    options                               = Helper.convertstruct2cell(OptionsStruct);
    
    solver                                = VariationalSolver2D.DipolarGas(options{:});
    pot                                   = VariationalSolver2D.Potentials(options{:});
    solver.Potential                      = pot.trap();

    % Run Solver
    [Params, Transf, psi, V, VDk]         = solver.run();    
    
end