From 54fc166eb419540c77d86bef7302b1d92447b267 Mon Sep 17 00:00:00 2001
From: Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
Date: Thu, 9 Jul 2020 10:56:59 +0200
Subject: [PATCH] Small changes in the SPRKKR STO.py example.

---
 tests/sprkkr/STO/STO.py | 73 ++++++++++++++++++++++++-----------------
 1 file changed, 43 insertions(+), 30 deletions(-)

diff --git a/tests/sprkkr/STO/STO.py b/tests/sprkkr/STO/STO.py
index 75ebb73..6bd5a18 100644
--- a/tests/sprkkr/STO/STO.py
+++ b/tests/sprkkr/STO/STO.py
@@ -11,6 +11,7 @@ import numpy as np
 import logging
 import sys
 import os
+import shutil
 import glob
 
 
@@ -78,8 +79,8 @@ if 'msspec' in sys.argv:
     # Build the SrTiO3 cluster                                                #
     ###########################################################################
     # We start by loading the potential since it updates info in the STO cell
-    pot = SPRKKRPotential(STO, "SrTiO3/SrTiO3.pot_new",
-                          *glob.glob("SrTiO3/*PHAGEN.pot"))
+    pot = SPRKKRPotential(STO, f"{prefix}/{prefix}.pot_new",
+                          *glob.glob(f"{prefix}/*PHAGEN.pot"))
 
     # tag each atom in the STO object to easily set the emitter in the
     # hemispherical_cluster function
@@ -94,7 +95,7 @@ if 'msspec' in sys.argv:
                                     emitter_plane=1, emitter_tag=1)
 
     # The created cluster is centered on the emitter atom, so defining
-    # the absorber attribute is straight forward:
+    # the absorber attribute is straightforward:
     cluster.absorber = get_atom_index(cluster, 0, 0, 0)
 
     ###########################################################################
@@ -105,37 +106,49 @@ if 'msspec' in sys.argv:
 
     # minimalistic set of parameter for XPD scan
     calc.set_atoms(cluster)
-    calc.source_parameters.energy = 700
     calc.calculation_parameters.scattering_order = 2
+    # We need to impose a maximum number of tl to use because there is still
+    # a memory bug that I need to investigate. Anyway, usually 25 is not that
+    # bad for the kind of atoms and the photon energy of lab sources...
+    calc.tmatrix_parameters.lmax_mode = 'imposed'
+    calc.tmatrix_parameters.lmaxt = 25
 
-    # Run a polar scan with the default MuffinTin potential for Ti(2p3/2)
-    data = calc.get_theta_scan(level='2p3/2')
+    data = None
+    for source_energy in np.arange(500., 1501., 100.):
+        calc.source_parameters.energy = source_energy
+        # Run a polar scan with the default MuffinTin potential for Ti(2p3/2)
+        calc.tmatrix_parameters.potential = 'muffin_tin'
+        data = calc.get_theta_scan(level='2p3/2', data=data)
 
-    # Now we use the previously generated SPRKKR potential
-    # and run the same calculation with appending the data to the previous one
-    calc.tmatrix_parameters.potential = pot
-    data = calc.get_theta_scan(level='2p3/2', data=data)
+        # Now we use the previously generated SPRKKR potential and run the same
+        # calculation
+        calc.tmatrix_parameters.potential = pot
+        data = calc.get_theta_scan(level='2p3/2', data=data)
+
+        # To better see the differences, plot the normalized signal on the
+        # same graph
+        # add a new dataset for storing normalized values
+        dset = data.add_dset(f"comparison at {source_energy} eV")
+        # make a copy of previous scans values
+        theta, muffintin_cs, sprkkr_cs = (data[-3].theta.copy(),
+                                        data[-3].cross_section.copy(),
+                                        data[-2].cross_section.copy())
+        # divide by the max
+        sprkkr_cs /= sprkkr_cs.max()
+        muffintin_cs /= muffintin_cs.max()
+        # add a new dataset with those values
+        dset.add_columns(theta=theta, sprkkr=sprkkr_cs, muffintin=muffintin_cs)
+        # add a view for this dataset
+        view = dset.add_view('Comparison',
+                            title=(f'Comparison of XPD polar scans for '
+                                   r'Ti(2p3/2) at $h\nu$ = {:.0f} eV'.format(
+                                       source_energy)),
+                            xlabel=r'$\Theta (\degree$)',
+                            ylabel='Normalized Signal (a.u.)')
+        view.select('theta', 'sprkkr', legend='With SPRKKR potential')
+        view.select('theta', 'muffintin', legend='With internal MT potential')
+        view.set_plot_options(autoscale=True)
 
-    # To better see the differences, plot the normalized signal on the
-    # same graph
-    # pick up previous data
-    theta, muffintin_cs, sprkkr_cs = (data[0].theta.copy(),
-                                      data[0].cross_section.copy(),
-                                      data[1].cross_section.copy())
-    # divide by the max
-    sprkkr_cs /= sprkkr_cs.max()
-    muffintin_cs /= muffintin_cs.max()
-    # add a new dataset with those values
-    dset = data.add_dset('comparison')
-    dset.add_columns(theta=theta, sprkkr=sprkkr_cs, muffintin=muffintin_cs)
-    # add a view for this dataset
-    view = dset.add_view('Comparison',
-                         title=r'Comparison of XPD polar scans for Ti(2p3/2)',
-                         xlabel=r'$\Theta (\degree$)',
-                         ylabel='Normalized Signal (a.u.)')
-    view.select('theta', 'sprkkr', legend='With SPRKKR potential')
-    view.select('theta', 'muffintin', legend='With internal MuffinTin potential')
-    view.set_plot_options(autoscale=True)
 
     # Pop up the final result
     data.view()