epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Add LEED calculator 

The Low Energy Electron Diffraction spectroscopy is now included
with Matrix Inversion algoriyhm only for now.
This commit is contained in:
Sylvain Tricot 2024-01-25 09:48:44 +01:00
parent 6c7038cdde
commit 7b104d5b20
17 changed files with 12544 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

413
src/msspec/calculator.py
View File

@ -17,8 +17,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/calculator.py
# Last modified: Tue, 25 Oct 2022 16:21:38 +0200
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1666707698 +0200
# Last modified: Thu, 25 Jan 2024 09:48:44 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
"""
@ -83,6 +83,7 @@ from msspec.parameters import CompCurveGeneralParameters
from msspec.parameters import DetectorParameters
from msspec.parameters import EIGParameters
from msspec.parameters import GlobalParameters
from msspec.parameters import LEEDParameters
from msspec.parameters import MuffintinParameters
from msspec.parameters import PEDParameters
from msspec.parameters import PhagenMallocParameters
@ -98,6 +99,7 @@ from msspec.spec.fortran import _eig_pw
from msspec.spec.fortran import _phd_mi_noso_nosp_nosym
from msspec.spec.fortran import _phd_se_noso_nosp_nosym
from msspec.spec.fortran import _phd_ce_noso_nosp_nosym
from msspec.spec.fortran import _led_mi_noso_nosp_nosym
from msspec.spec.fortran import _comp_curves
from msspec.utils import get_atom_index
@ -150,6 +152,9 @@ class _MSCALCULATOR(Calculator):
if spectroscopy == 'PED':
self.spectroscopy_parameters = PEDParameters(self.phagen_parameters,
self.spec_parameters)
elif spectroscopy == 'LEED':
self.spectroscopy_parameters = LEEDParameters(self.phagen_parameters,
self.spec_parameters)
elif spectroscopy == 'EIG':
self.spectroscopy_parameters = EIGParameters(self.phagen_parameters,
self.spec_parameters)
@ -413,6 +418,14 @@ class _MSCALCULATOR(Calculator):
"an allowed combination.".format(self.global_parameters.spectroscopy,
self.global_parameters.algorithm))
raise ValueError
elif self.global_parameters.spectroscopy == 'LEED':
if self.global_parameters.algorithm == 'inversion':
do_spec = _led_mi_noso_nosp_nosym.run
else:
LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
"an allowed combination.".format(self.global_parameters.spectroscopy,
self.global_parameters.algorithm))
raise ValueError
elif self.global_parameters.spectroscopy == 'EIG':
if self.global_parameters.algorithm == 'inversion':
do_spec = _eig_mi.run
@ -1047,6 +1060,400 @@ class _EIG(_MSCALCULATOR):
return self.iodata
class _LEED(_MSCALCULATOR):
"""This class creates a calculator object for Low Electron Energy Diffraction
spectroscopy.
:param algorithm: The algorithm to use for the computation. See
:ref:`globalparameters-algorithm` for more details about the allowed
values and the type.
:param polarization: The incoming light polarization (see
:ref:`globalparameters-polarization`)
:param dichroism: Wether to enable or not the dichroism (see
:ref:`globalparameters-dichroism`)
:param spinpol: Enable or disable the spin polarization in the calculation
(see :ref:`globalparameters-spinpol`)
:param folder: The path to the temporary folder for the calculations. See
:ref:`globalparameters-folder`
:param txt: The name of a file where to redirect standard output. The string
'-' will redirect the standard output to the screen (default).
:type txt: str
.. note::
This class constructor is not meant to be called directly by the user.
Use the :py:func:`MSSPEC` to instanciate any calculator.
"""
def __init__(self, algorithm='inversion', polarization=None, dichroism=None,
spinpol=False, folder='./calc', txt='-'):
_MSCALCULATOR.__init__(self, spectroscopy='LEED', algorithm=algorithm,
polarization=polarization, dichroism=dichroism,
spinpol=spinpol, folder=folder, txt=txt)
self.source_parameters.theta = 0
self.source_parameters.phi = 0
self.iodata = iodata.Data('LEED Simulation')
def _get_scan(self, scan_type='theta', phi=0,
theta=np.linspace(-70, 70, 141),
kinetic_energy=None, data=None,
malloc={}, other_parameters={}):
LOGGER.info("Computting the %s scan...", scan_type)
# Force absorber to be 0.
self.atoms.absorber = get_atom_index(self.atoms, 0, 0, 0)
self.detector_parameters.rotate = True
self.source_parameters.theta = 0
self.source_parameters.phi = 0
if data:
self.iodata = data
if kinetic_energy is None:
LOGGER.error('The kinetic energy is not specified!')
raise ValueError('You must define a kinetic_energy value.')
# update the parameters
self.scan_parameters.set_parameter('kinetic_energy', kinetic_energy)
all_ke = self.scan_parameters.get_parameter('ke_array')
if np.any(all_ke.value < 0):
LOGGER.error('Source energy is not high enough or level too deep!')
raise ValueError('Kinetic energy is < 0! ({})'.format(
kinetic_energy))
self.scan_parameters.set_parameter('type', scan_type)
# make sure there is only one energy point in scatf scan
if scan_type == 'scatf':
assert len(all_ke) == 1, ('kinetic_energy should not be an array '
'in scatf scan')
if scan_type != 'scatf':
self.scan_parameters.set_parameter('phi', phi)
self.scan_parameters.set_parameter('theta', theta)
#self.spectroscopy_parameters.set_parameter('level', level)
# It is still possible to modify any option right before runing phagen
# and spec
for k, v in other_parameters.items():
grp_str, param_str = k.split('.')
grp = getattr(self, grp_str)
grp.set_parameter(param_str, v, force=True)
self.get_tmatrix()
self.run_spec(malloc)
# Now load the data
ndset = len(self.iodata)
dset = self.iodata.add_dset('{} scan [{:d}]'.format(scan_type, ndset))
for p in self.get_parameters():
bundle = {'group': str(p.group),
'name': str(p.name),
'value': str(p.value),
'unit': '' if p.unit is None else str(p.unit)}
dset.add_parameter(**bundle)
if scan_type in ('theta', 'phi', 'energy'):
results_fname = os.path.join(self.tmp_folder, 'output/results.dat')
data = self.specio.load_results(results_fname)
for _plane, _theta, _phi, _energy, _dirsig, _cs in data.T:
if _plane == -1:
dset.add_row(theta=_theta, phi=_phi, energy=_energy, cross_section=_cs, direct_signal=_dirsig)
elif scan_type in ('scatf',):
results_fname = os.path.join(self.tmp_folder, 'output/facdif1.dat')
data = self.specio.load_facdif(results_fname)
data = data[:, [1, 4, 5, 6, 8]].T
_proto, _sf_real, _sf_imag, _theta, _energy = data
_sf = _sf_real + _sf_imag * 1j
dset.add_columns(proto_index=_proto, sf_real=np.real(_sf),
sf_imag=np.imag(_sf), sf_module=np.abs(_sf),
theta=_theta, energy=_energy)
elif scan_type in ('theta_phi',):
results_fname = os.path.join(self.tmp_folder, 'output/results.dat')
data = self.specio.load_results(results_fname)
#theta_c, phi_c = data[[2, 3], :]
#xsec_c = data[-1, :]
#dirsig_c = data[-2, :]
#dset.add_columns(theta=theta_c)
#dset.add_columns(phi=phi_c)
#dset.add_columns(cross_section=xsec_c)
#dset.add_columns(direct_signal=dirsig_c)
for _plane, _theta, _phi, _energy, _dirsig, _cs in data.T:
if _plane == -1:
dset.add_row(theta=_theta, phi=_phi, energy=_energy, cross_section=_cs,
direct_signal=_dirsig)
# create a view
title = ''
for ke in all_ke.value:
if scan_type == 'theta':
title = 'Polar scan at {:.2f} eV'.format(ke)
xlabel = r'Angle $\theta$($\degree$)'
ylabel = r'Signal (a. u.)'
view = dset.add_view("E = {:.2f} eV".format(ke), title=title,
xlabel=xlabel, ylabel=ylabel, autoscale=True)
for angle_phi in self.scan_parameters.get_parameter(
'phi').value:
where = ("energy=={:.2f} and phi=={:.2f}"
"").format(ke, angle_phi)
legend = r'$\phi$ = {:.1f} $\degree$'.format(angle_phi)
view.select('theta', 'cross_section', where=where,
legend=legend)
if scan_type == 'phi':
title = 'Azimuthal scan at {:.2f} eV'.format(ke)
xlabel = r'Angle $\phi$($\degree$)'
ylabel = r'Signal (a. u.)'
view = dset.add_view("E = {:.2f} eV".format(ke), title=title,
xlabel=xlabel, ylabel=ylabel)
for angle_theta in self.scan_parameters.get_parameter(
'theta').value:
where = ("energy=={:.2f} and theta=={:.2f}"
"").format(ke, angle_theta)
legend = r'$\theta$ = {:.1f} $\degree$'.format(angle_theta)
view.select('phi', 'cross_section', where=where,
legend=legend)
if scan_type == 'theta_phi':
absorber_symbol = self.atoms[self.atoms.absorber].symbol
title = ('Stereographic projection at {:.2f} eV'
'').format(ke)
xlabel = r'Angle $\phi$($\degree$)'
ylabel = r'Signal (a. u.)'
view = dset.add_view("E = {:.2f} eV".format(ke), title=title,
xlabel=xlabel, ylabel=ylabel,
projection='stereo', colorbar=True, autoscale=True)
view.select('theta', 'phi', 'cross_section')
if scan_type == 'scatf':
for i in range(self.phagenio.nat):
proto_index = i+1
title = 'Scattering factor at {:.3f} eV'.format(kinetic_energy)
mini = min(map(np.min, [dset.sf_real, dset.sf_imag, dset.sf_module]))
maxi = max(map(np.max, [dset.sf_real, dset.sf_imag, dset.sf_module]))
view = dset.add_view("Proto. atom #{:d}".format(proto_index),
title=title, projection='polar',
ylim=[mini, maxi])
where = "proto_index=={:d}".format(proto_index)
view.select('theta', 'sf_module', where=where,
legend=r'$|f(\theta)|$')
view.select('theta', 'sf_real', where=where,
legend=r'$\Re(f(\theta))$')
view.select('theta', 'sf_imag', where=where,
legend=r'$\Im(f(\theta))$')
if scan_type == 'energy':
absorber_symbol = self.atoms[self.atoms.absorber].symbol
title = (r'Energy scan of {}({}) at $\theta$={:.2f}$\degree$ and '
'$\phi$={:.2f}$\degree$').format(
absorber_symbol, level, theta, phi)
xlabel = r'Photoelectron kinetic energy (eV)'
ylabel = r'Signal (a. u.)'
view = dset.add_view("EnergyScan".format(ke), title=title,
xlabel=xlabel, ylabel=ylabel)
view.select('energy', 'cross_section')
# save the cluster
#clusbuf = StringIO()
#self.atoms.info['absorber'] = self.atoms.absorber
#self.atoms.write(clusbuf, format='xyz')
#dset.add_parameter(group='Cluster', name='cluster', value=clusbuf.getvalue(), hidden="True")
self.add_cluster_to_dset(dset)
LOGGER.info('%s scan computing done!', scan_type)
return self.iodata
def get_potential(self, atom_index=None, data=None, units={'energy': 'eV', 'space': 'angstrom'}):
"""Computes the coulombic part of the atomic potential.
:param atom_index: The atom indices to get the potential of, either as a list or as a single integer
:param data: The data object to store the results to
:param units: The units to be used. A dictionary with the keys 'energy' and 'space'
:return: A Data object
"""
LOGGER.info("Getting the Potential...")
LOGGER.debug(get_call_info(inspect.currentframe()))
_units = {'energy': 'eV', 'space': 'angstrom'}
_units.update(units)
if data:
self.iodata = data
self.run_phagen()
filename = os.path.join(self.tmp_folder, 'output/tmatrix.tl')
tl = self.phagenio.load_tl_file(filename)
filename = os.path.join(self.tmp_folder, 'output/cluster.clu')
self.phagenio.load_cluster_file(filename)
if self.phagen_parameters.potgen in ('in'):
filename = os.path.join(self.tmp_folder, 'output/plot/plot_vc.dat')
else:
filename = os.path.join(self.tmp_folder, 'output/plot/plot_v.dat')
pot_data = self.phagenio.load_potential_file(filename)
cluster = self.phagen_parameters.get_parameter('atoms').value
dset = self.iodata.add_dset('Potential [{:d}]'.format(len(self.iodata)))
r = []
v = []
index = np.empty((0,1), dtype=int)
absorber_position = cluster[cluster.absorber].position
for _pot_data in pot_data:
# find the proto index of these data
at_position = (_pot_data['coord'] * UREG.bohr_radius).to('angstrom').magnitude + absorber_position
at_index = get_atom_index(cluster, *at_position)
at_proto_index = cluster[at_index].get('proto_index')
#values = np.asarray(_pot_data['values'])
values = _pot_data['values']
index = np.append(index, np.ones(values.shape[0], dtype=int) * at_proto_index)
r = np.append(r, (values[:, 0] * UREG.bohr_radius).to(_units['space']).magnitude)
v = np.append(v, (values[:, 1] * UREG.rydberg).to(_units['energy']).magnitude)
dset.add_columns(distance=r, potential=v, index=index)
view = dset.add_view('potential data', title='Potential energy of atoms',
xlabel='distance from atomic center [{:s}]'.format(_units['space']),
ylabel='energy [{:s}]'.format(_units['energy']), scale='linear',
autoscale=True)
if atom_index == None:
for i in range(pot_data[len(pot_data) - 1]['index']):
view.select('distance', 'potential', where="index=={:d}".format(i),
legend="Atom index #{:d}".format(i + 1))
else:
for i in atom_index:
view.select('distance', 'potential', where="index=={:d}".format(cluster[i].get('proto_index') - 1),
legend="Atom index #{:d}".format(i))
return self.iodata
def get_scattering_factors(self, level='1s', kinetic_energy=None,
data=None, **kwargs):
"""Computes the scattering factors of all prototypical atoms in the
cluster.
This function computes the real and imaginery parts of the scattering
factor as well as its modulus for each non symetrically equivalent atom
in the cluster. The results are stored in the *data* object if provided
as a parameter.
:param level: The electronic level. See :ref:`pedparameters-level`.
:param kinetic_energy: see :ref:`scanparameters-kinetic_energy`.
:param data: a :py:class:`iodata.Data` object to append the results to
or None.
:returns: The modified :py:class:`iodata.Data` object passed as an
argument or a new :py:class:`iodata.Data` object.
"""
data = self._get_scan(scan_type='scatf', level=level, data=data,
kinetic_energy=kinetic_energy, **kwargs)
return data
def get_theta_scan(self, phi=0, theta=np.linspace(-70, 70, 141),
kinetic_energy=None, data=None, **kwargs):
"""Computes a polar scan of the emitted photoelectrons.
:param phi: The azimuthal angle in degrees. See
:ref:`scanparameters-phi`.
:param theta: All the values of the polar angle to be computed. See
:ref:`scanparameters-theta`.
:param level: The electronic level. See :ref:`pedparameters-level`.
:param kinetic_energy: see :ref:`scanparameters-kinetic_energy`.
:param data: a :py:class:`iodata.Data` object to append the results to
or None.
:returns: The modified :py:class:`iodata.Data` object passed as an
argument or a new :py:class:`iodata.Data` object.
"""
data = self._get_scan(scan_type='theta', theta=theta,
phi=phi, kinetic_energy=kinetic_energy,
data=data, **kwargs)
return data
def get_phi_scan(self, phi=np.linspace(0, 359, 359), theta=0,
kinetic_energy=None, data=None, **kwargs):
"""Computes an azimuthal scan of the emitted photoelectrons.
:param phi: All the values of the azimuthal angle to be computed. See
:ref:`scanparameters-phi`.
:param theta: The polar angle in degrees. See
:ref:`scanparameters-theta`.
:param level: The electronic level. See :ref:`pedparameters-level`.
:param kinetic_energy: see :ref:`scanparameters-kinetic_energy`.
:param data: a :py:class:`iodata.Data` object to append the results to
or None.
:returns: The modified :py:class:`iodata.Data` object passed as an
argument or a new :py:class:`iodata.Data` object.
"""
data = self._get_scan(scan_type='phi', theta=theta,
phi=phi, kinetic_energy=kinetic_energy,
data=data, **kwargs)
return data
def get_theta_phi_scan(self, phi=np.linspace(0, 360),
theta=np.linspace(0, 90, 45),
kinetic_energy=None, data=None, **kwargs):
"""Computes a stereographic scan of the emitted photoelectrons.
The azimuth ranges from 0 to 360° and the polar angle ranges from 0 to
90°.
:param level: The electronic level. See :ref:`pedparameters-level`.
:param kinetic_energy: see :ref:`scanparameters-kinetic_energy`.
:param data: a :py:class:`iodata.Data` object to append the results to
or None.
:returns: The modified :py:class:`iodata.Data` object passed as an
argument or a new :py:class:`iodata.Data` object.
"""
data = self._get_scan(scan_type='theta_phi', theta=theta,
phi=phi, kinetic_energy=kinetic_energy, data=data,
**kwargs)
return data
def get_energy_scan(self, phi=0, theta=0,
level=None, kinetic_energy=None, data=None, **kwargs):
"""Computes an energy scan of the emitted photoelectrons.
:param phi: All the values of the azimuthal angle to be computed. See
:ref:`scanparameters-phi`.
:param theta: The polar angle in degrees. See
:ref:`scanparameters-theta`.
:param level: The electronic level. See :ref:`pedparameters-level`.
:param kinetic_energy: see :ref:`scanparameters-kinetic_energy`.
:param data: a :py:class:`iodata.Data` object to append the results to
or None.
:returns: The modified :py:class:`iodata.Data` object passed as an
argument or a new :py:class:`iodata.Data` object.
"""
data = self._get_scan(scan_type='energy', level=level, theta=theta,
phi=phi, kinetic_energy=kinetic_energy,
data=data, **kwargs)
return data
def MSSPEC(spectroscopy='PED', **kwargs):
""" The MsSpec calculator constructor.
@ -1285,5 +1692,7 @@ class RFACTOR(object):
if __name__ == "__main__":
pass

37
src/msspec/iodata.py
View File

@ -17,8 +17,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/iodata.py
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
# Last modified: Thu, 25 Jan 2024 09:48:44 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
"""
@ -86,6 +86,7 @@ from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.backends.backend_wxagg import NavigationToolbar2WxAgg
from matplotlib.figure import Figure
from matplotlib.ticker import FormatStrFormatter
from terminaltables import AsciiTable
import msspec
@ -93,13 +94,17 @@ from msspec.msspecgui.msspec.gui.clusterviewer import ClusterViewer
from msspec.misc import LOGGER
def cols2matrix(x, y, z, nx=88*1+1, ny=360*1+1):
def cols2matrix(x, y, z, nx=88*1+1, ny=360*1+1, xlim=[None, None], ylim=[None, None]):
# mix the values of existing theta and new theta and return the
# unique values
newx = np.linspace(np.min(x), np.max(x), nx)
xmin = xlim[0] if xlim[0] is not None else np.min(x)
xmax = xlim[1] if xlim[1] is not None else np.max(x)
ymin = ylim[0] if ylim[0] is not None else np.min(y)
ymax = ylim[1] if ylim[1] is not None else np.max(y)
newx = np.linspace(xmin, xmax, nx)
newy = np.linspace(np.min(y), np.max(y), ny)
ux = np.unique(np.append(x, newx))
uy = np.unique(np.append(y, newy))
ux = np.unique(np.sort(np.append(x, newx)).clip(xmin, xmax))
uy = np.unique(np.sort(np.append(y, newy)).clip(ymin, ymax))
# create an empty matrix to hold the results
zz = np.empty((len(ux), len(uy)))
@ -813,7 +818,8 @@ class _DataSetView(object):
xlabel='', ylabel='', grid=True, legend=[], colorbar=False,
projection='rectilinear', xlim=[None, None], ylim=[None, None],
scale='linear',
marker=None, autoscale=False)
specular=None,
marker=None, autoscale=True)
self._plotopts.update(plotopts)
self._selection_tags = []
self._selection_conditions = []
@ -879,19 +885,26 @@ class _DataSetView(object):
axes.set_xticks(xvalues)
else:
if proj in ('ortho', 'stereo'):
theta, phi, Xsec = cols2matrix(*values)
theta_ticks = np.arange(0, 91, 15)
theta, phi, Xsec = cols2matrix(*values, xlim=opts['xlim'], ylim=opts['ylim'])
#theta_ticks = np.arange(0, 91, 15)
theta_ticks = np.linspace(np.min(theta), np.max(theta), 7)
if proj == 'ortho':
R = np.sin(np.radians(theta))
R_ticks = np.sin(np.radians(theta_ticks))
elif proj == 'stereo':
R = 2 * np.tan(np.radians(theta/2.))
R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
#R = np.tan(np.radians(theta/2.))
X, Y = np.meshgrid(np.radians(phi), R)
if opts['specular'] is not None:
Xsec[Y<np.radians(opts['specular'])] = np.nan
im = axes.pcolormesh(X, Y, Xsec, shading='gouraud')
axes.set_yticks(R_ticks)
axes.set_yticklabels(theta_ticks)
axes.set_yticklabels([ '{:.1f}'.format(_) for _ in theta_ticks ])
#print(R_ticks)
print(theta_ticks)
#exit()
figure.colorbar(im)
@ -914,8 +927,6 @@ class _DataSetView(object):
axes.set_title(opts['title'])
axes.set_xlabel(opts['xlabel'])
axes.set_ylabel(opts['ylabel'])
axes.set_xlim(*opts['xlim'])
axes.set_ylim(*opts['ylim'])
#axes.set_pickradius(5)
if label:
axes.legend()

69
src/msspec/parameters.py
View File

@ -19,8 +19,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/parameters.py
# Last modified: Tue, 15 Feb 2022 15:37:28 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
# Last modified: Thu, 25 Jan 2024 09:48:45 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
"""
@ -488,11 +488,11 @@ class SpecParameters(BaseParameters):
fmt='.2f'),
Parameter('leed_r1', types=float, default=-1.0,
fmt='.3f'),
Parameter('leed_thini', types=float, default=-55.0,
Parameter('leed_thini', types=float, default=0.,
fmt='.2f'),
Parameter('leed_phiini', types=float, default=0.,
fmt='.2f'),
Parameter('leed_imod', types=int, default=1,
Parameter('leed_imod', types=int, default=0,
fmt='d'),
Parameter('leed_imoy', types=int, default=0,
fmt='d'),
@ -833,7 +833,7 @@ class GlobalParameters(BaseParameters):
'AED': ('aed', 'AED'),
'LEED': ('led', 'LED'),
'EXAFS': ('xas', 'XAS'),
'EIG': ('xpd', 'EIG'),
'EIG': ('led', 'EIG'),
}
phagen_calctype, spec_calctype = mapping[p.value]
self.phagen_parameters.calctype = phagen_calctype
@ -1154,7 +1154,7 @@ class DetectorParameters(BaseParameters):
default=None, doc="""
Used to averaged the signal over directions lying in the
cone of half-angle *angular_acceptance*. The number of
directions to take into account depends on the choosen
directions to take into account depends on the chosen
value:
- **None**, for no averaging at all
@ -1307,7 +1307,7 @@ class ScanParameters(BaseParameters):
calculation_parameters.set_parameter('basis_functions',
'spherical', force=True)
LOGGER.info('\'%s\' scan type choosen.', p.value)
LOGGER.info('\'%s\' scan type chosen.', p.value)
def bind_theta(self, p):
spectro = self.global_parameters.spectroscopy
@ -1422,6 +1422,10 @@ class ScanParameters(BaseParameters):
self.spec_parameters.ped_e1 = energies[1]
self.spec_parameters.ped_ne = energies[2]
self.spec_parameters.leed_e0 = energies[0]
self.spec_parameters.leed_e1 = energies[1]
self.spec_parameters.leed_ne = energies[2]
self.spec_parameters.eigval_ekini = energies[0]
self.spec_parameters.eigval_ekfin = energies[1]
self.spec_parameters.eigval_ne = energies[2]
@ -1829,6 +1833,57 @@ class PEDParameters(BaseParameters):
self.spec_parameters.ped_iso = somap[p.value]
class LEEDParameters(BaseParameters):
def __init__(self, phagen_parameters, spec_parameters):
# parameters = (
# Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$',
# default='1s', doc="""
# The level is the electronic level where the electron comes from.
# It is written: *nlJ*
# where:
# - *n* is the principal quantum number
# - *l* is the orbital quantum number
# - *J* is the spin-orbit component
# Example::
# >>> calc.spectroscopy_parameters.level = '2p3/2'
# >>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2'
# """),
# Parameter('final_state', types=int, limits=(-1, 2), default=2),
# Parameter('spin_orbit', types=(type(None), str),
# allowed_values=(None, 'single', 'both'), default=None),
# )
BaseParameters.__init__(self)
#self.add_parameters(*parameters)
self.phagen_parameters = phagen_parameters
self.spec_parameters = spec_parameters
# def bind_level(self, p):
# edge = get_level_from_electron_configuration(p.value)
# self.phagen_parameters.edge = edge
# li, so = re.match(r'(^\d+[spdfg])(.*$)', p.value).groups()
# if so == '':
# so = '1/2'
# self.spec_parameters.ped_li = li
# self.spec_parameters.ped_so = so
# self.spec_parameters.extra_level = p.value
# def bind_final_state(self, p):
# self.spec_parameters.ped_initl = p.value
# def bind_spin_orbit(self, p):
# somap = {
# None: 0,
# 'single': 1,
# 'both': 2}
# self.spec_parameters.ped_iso = somap[p.value]
class EIGParameters(BaseParameters):
def __init__(self, phagen_parameters, spec_parameters):
parameters = (

8
src/msspec/spec/fortran/Makefile
View File

@ -1,6 +1,6 @@
.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw comp_curve clean
.PHONY: all phd_se phd_mi phd_ce led_mi eig_mi eig_pw comp_curve clean
all: phd_se phd_mi phd_ce eig_mi eig_pw comp_curve
all: phd_se phd_mi phd_ce led_mi eig_mi eig_pw comp_curve
phd_se:
@+$(MAKE) -f phd_se_noso_nosp_nosym.mk all
@ -11,6 +11,9 @@ phd_mi:
phd_ce:
@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk all
led_mi:
@+$(MAKE) -f led_mi_noso_nosp_nosym.mk all
eig_mi:
@+$(MAKE) -f eig_mi.mk all
@ -24,6 +27,7 @@ clean::
@+$(MAKE) -f phd_se_noso_nosp_nosym.mk $@
@+$(MAKE) -f phd_mi_noso_nosp_nosym.mk $@
@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk $@
@+$(MAKE) -f led_mi_noso_nosp_nosym.mk $@
@+$(MAKE) -f eig_mi.mk $@
@+$(MAKE) -f eig_pw.mk $@
@+$(MAKE) -f comp_curve.mk $@

11
src/msspec/spec/fortran/led_mi_noso_nosp_nosym.mk Normal file
View File

@ -0,0 +1,11 @@
memalloc_src := memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f
cluster_gen_src := $(wildcard cluster_gen/*.f)
common_sub_src := $(wildcard common_sub/*.f)
renormalization_src := $(wildcard renormalization/*.f)
led_mi_noso_nosp_nosym_src := $(filter-out led_mi_noso_nosp_nosym/lapack_inv.f, $(wildcard led_mi_noso_nosp_nosym/*.f))
SRCS = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(led_mi_noso_nosp_nosym_src)
MAIN_F = led_mi_noso_nosp_nosym/main.f
SO = _led_mi_noso_nosp_nosym.so
include ../../../options.mk

85
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/dwsph.f Normal file
View File

@ -0,0 +1,85 @@
C
C=======================================================================
C
SUBROUTINE DWSPH(JTYP,JE,X,TLT,ISPEED)
C
C This routine recomputes the T-matrix elements taking into account the
C mean square displacements.
C
C When the argument X is tiny, no vibrations are taken into account
C
C Last modified : 25 Apr 2013
C
USE DIM_MOD
C
USE TRANS_MOD
C
DIMENSION GNT(0:N_GAUNT)
C
COMPLEX TLT(0:NT_M,4,NATM,NE_M),SL1,ZEROC
C
COMPLEX*16 FFL(0:2*NL_M)
C
DATA PI4,EPS /12.566371,1.0E-10/
C
ZEROC=(0.,0.)
C
IF(X.GT.EPS) THEN
C
C Standard case: vibrations
C
IF(ISPEED.LT.0) THEN
NSUM_LB=ABS(ISPEED)
ENDIF
C
COEF=PI4*EXP(-X)
NL2=2*LMAX(JTYP,JE)+2
IBESP=5
MG1=0
MG2=0
C
CALL BESPHE(NL2,IBESP,X,FFL)
C
DO L=0,LMAX(JTYP,JE)
XL=FLOAT(L+L+1)
SL1=ZEROC
C
DO L1=0,LMAX(JTYP,JE)
XL1=FLOAT(L1+L1+1)
CALL GAUNT(L,MG1,L1,MG2,GNT)
L2MIN=ABS(L1-L)
IF(ISPEED.GE.0) THEN
L2MAX=L1+L
ELSEIF(ISPEED.LT.0) THEN
L2MAX=L2MIN+2*(NSUM_LB-1)
ENDIF
SL2=0.
C
DO L2=L2MIN,L2MAX,2
XL2=FLOAT(L2+L2+1)
C=SQRT(XL1*XL2/(PI4*XL))
SL2=SL2+C*GNT(L2)*REAL(DREAL(FFL(L2)))
ENDDO
C
SL1=SL1+SL2*TL(L1,1,JTYP,JE)
ENDDO
C
TLT(L,1,JTYP,JE)=COEF*SL1
C
ENDDO
C
ELSE
C
C Argument X tiny: no vibrations
C
DO L=0,LMAX(JTYP,JE)
C
TLT(L,1,JTYP,JE)=TL(L,1,JTYP,JE)
C
ENDDO
C
ENDIF
C
RETURN
C
END

26
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/facdif.f Normal file
View File

@ -0,0 +1,26 @@
C
C=======================================================================
C
SUBROUTINE FACDIF(COSTH,JAT,JE,FTHETA)
C
C This routine computes the plane wave scattering factor
C
USE DIM_MOD
C
USE TRANS_MOD
C
DIMENSION PL(0:100)
C
COMPLEX FTHETA
C
FTHETA=(0.,0.)
NL=LMAX(JAT,JE)+1
CALL POLLEG(NL,COSTH,PL)
DO 20 L=0,NL-1
FTHETA=FTHETA+(2*L+1)*TL(L,1,JAT,JE)*PL(L)
20 CONTINUE
FTHETA=FTHETA/VK(JE)
C
RETURN
C
END

113
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/facdif1.f Normal file
View File

@ -0,0 +1,113 @@
C
C=======================================================================
C
SUBROUTINE FACDIF1(VKE,RJ,RJK,THRJ,PHIRJ,BETA,GAMMA,L,M,FSPH,JAT,J
&E,*)
C
C This routine computes a spherical wave scattering factor
C
C Last modified : 03/04/2006
C
USE DIM_MOD
USE APPROX_MOD
USE EXPFAC_MOD
USE TRANS_MOD
USE TYPCAL_MOD , I2 => IPHI, I3 => IE, I4 => ITHETA, I5 => IMOD, I
&6 => IPOL, I7 => I_CP, I8 => I_EXT, I9 => I_TEST
C
DIMENSION PLMM(0:100,0:100)
DIMENSION D(1-NL_M:NL_M-1,1-NL_M:NL_M-1,0:NL_M-1)
C
COMPLEX HLM(0:NO_ST_M,0:NL_M-1),HLN(0:NO_ST_M,0:NL_M-1),FSPH,RHOJ
COMPLEX HLM1,HLM2,HLM3,HLM4,ALMU,BLMU,SLP,SNU,SMU,VKE
COMPLEX RHOJK
C
C
DATA PI/3.141593/
C
A=1.
INTER=0
IF(ITL.EQ.1) VKE=VK(JE)
RHOJ=VKE*RJ
RHOJK=VKE*RJK
HLM1=(1.,0.)
HLM2=(1.,0.)
HLM3=(1.,0.)
HLM4=(1.,0.)
IEM=1
CSTH=COS(BETA)
IF((IFTHET.EQ.0).OR.(THRJ.LT.0.0001)) THEN
INTER=1
BLMU=SQRT(4.*PI/FLOAT(2*L+1))*CEXP((0.,-1.)*M*(PHIRJ-PI))
ENDIF
CALL PLM(CSTH,PLMM,LMAX(JAT,JE))
IF(ISPHER.EQ.0) NO1=0
IF(ISPHER.EQ.1) THEN
IF(NO.EQ.8) THEN
NO1=LMAX(JAT,JE)+1
ELSE
NO1=NO
ENDIF
CALL POLHAN(ISPHER,NO1,LMAX(JAT,JE),RHOJ,HLM)
IF(IEM.EQ.0) THEN
HLM4=HLM(0,L)
ENDIF
IF(RJK.GT.0.0001) THEN
NDUM=0
CALL POLHAN(ISPHER,NDUM,LMAX(JAT,JE),RHOJK,HLN)
ENDIF
CALL DJMN(THRJ,D,L)
A1=ABS(D(0,M,L))
IF(((A1.LT.0.0001).AND.(IFTHET.EQ.1)).AND.(INTER.EQ.0)) RETURN 1
&
ENDIF
MUMAX=MIN0(L,NO1)
SMU=(0.,0.)
DO 10 MU=0,MUMAX
IF(MOD(MU,2).EQ.0) THEN
B=1.
ELSE
B=-1.
IF(SIN(BETA).LT.0.) THEN
A=-1.
ENDIF
ENDIF
IF(ISPHER.LE.1) THEN
ALMU=(1.,0.)
C=1.
ENDIF
IF(ISPHER.EQ.0) GOTO 40
IF(INTER.EQ.0) BLMU=CMPLX(D(M,0,L))
IF(MU.GT.0) THEN
C=B*FLOAT(L+L+1)/EXPF(MU,L)
ALMU=(D(M,MU,L)*CEXP((0.,-1.)*MU*GAMMA)+B*
* CEXP((0.,1.)*MU*GAMMA)*D(M,-MU,L))/BLMU
ELSE
C=1.
ALMU=CMPLX(D(M,0,L))/BLMU
ENDIF
40 SNU=(0.,0.)
NU1=INT(0.5*(NO1-MU)+0.0001)
NUMAX=MIN0(NU1,L-MU)
DO 20 NU=0,NUMAX
SLP=(0.,0.)
LPMIN=MAX0(MU,NU)
DO 30 LP=LPMIN,LMAX(JAT,JE)
IF(ISPHER.EQ.1) THEN
HLM1=HLM(NU,LP)
IF(RJK.GT.0.0001) HLM3=HLN(0,LP)
ENDIF
SLP=SLP+FLOAT(2*LP+1)*TL(LP,1,JAT,JE)*HLM1*PLMM(LP,MU)*HLM3
30 CONTINUE
IF(ISPHER.EQ.1) THEN
HLM2=HLM(MU+NU,L)
ENDIF
SNU=SNU+SLP*HLM2
20 CONTINUE
SMU=SMU+SNU*C*ALMU*A*B
10 CONTINUE
FSPH=SMU/(VKE*HLM4)
C
RETURN
C
END

192
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/inv_mat_ms2_la.f Normal file
View File

@ -0,0 +1,192 @@
C
C=======================================================================
C
SUBROUTINE INV_MAT_MS(JE,TAU)
C
C This subroutine stores the multiple scattering matrix and invert
C it to obtain the scattering path operator exactly.
C
C (Photoelectron case)
C
C Last modified : 24 Apr 2007
C
C INCLUDE 'spec.inc'
USE DIM_MOD
USE COOR_MOD
USE INIT_L_MOD
USE TRANS_MOD
C
C PARAMETER(NLTWO=2*NL_M)
C
COMPLEX*16 HL1(0:2*NL_M),SM(LINMAX*NATCLU_M,LINMAX*NATCLU_M)
COMPLEX*16 SUM_L,ONEC,IC,ZEROC,WORK(4*LINMAX*NATCLU_M)
COMPLEX*16 YLM(0:2*NL_M,-2*NL_M:2*NL_M),TLJ,TLK,EXPKJ
C
COMPLEX TAU(LINMAX*LINMAX*NATCLU_M*NATCLU_M)
C
C
REAL*8 PI,ATTKJ,GNT(0:N_GAUNT),XKJ,YKJ,ZKJ,RKJ,ZDKJ,KRKJ
C
INTEGER IPIV(LINMAX*NATCLU_M)
C
C
DATA PI /3.1415926535898D0/
C
ONEC=(1.D0,0.D0)
IC=(0.D0,1.D0)
ZEROC=(0.D0,0.D0)
IBESS=3
C
C Construction of the multiple scattering matrix MS = (I-GoT).
C Elements are stored using a linear index LINJ representing
C (J,LJ)
C
JLIN=0
DO JTYP=1,N_PROT
NBTYPJ=NATYP(JTYP)
LMJ=LMAX(JTYP,JE)
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
XJ=SYM_AT(1,JATL)
YJ=SYM_AT(2,JATL)
ZJ=SYM_AT(3,JATL)
C
DO LJ=0,LMJ
ILJ=LJ*LJ+LJ+1
TLJ=DCMPLX(TL(LJ,1,JTYP,JE))
DO MJ=-LJ,LJ
INDJ=ILJ+MJ
JLIN=JLIN+1
C
KLIN=0
DO KTYP=1,N_PROT
NBTYPK=NATYP(KTYP)
LMK=LMAX(KTYP,JE)
DO KNUM=1,NBTYPK
KATL=NCORR(KNUM,KTYP)
IF(KATL.NE.JATL) THEN
XKJ=DBLE(SYM_AT(1,KATL)-XJ)
YKJ=DBLE(SYM_AT(2,KATL)-YJ)
ZKJ=DBLE(SYM_AT(3,KATL)-ZJ)
RKJ=DSQRT(XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ)
KRKJ=DBLE(VK(JE))*RKJ
ATTKJ=DEXP(-DIMAG(DCMPLX(VK(JE)))*
& RKJ)
EXPKJ=(XKJ+IC*YKJ)/RKJ
ZDKJ=ZKJ/RKJ
CALL SPH_HAR2(2*NL_M,ZDKJ,EXPKJ,YLM,
& LMJ+LMK)
CALL BESPHE2(LMJ+LMK+1,IBESS,KRKJ,
& HL1)
ENDIF
C
DO LK=0,LMK
ILK=LK*LK+LK+1
L_MIN=ABS(LK-LJ)
L_MAX=LK+LJ
TLK=DCMPLX(TL(LK,1,KTYP,JE))
DO MK=-LK,LK
INDK=ILK+MK
KLIN=KLIN+1
SM(KLIN,JLIN)=ZEROC
SUM_L=ZEROC
IF(KATL.NE.JATL) THEN
CALL GAUNT2(LK,MK,LJ,MJ,GNT)
C
DO L=L_MIN,L_MAX,2
M=MJ-MK
IF(ABS(M).LE.L) THEN
SUM_L=SUM_L+(IC**L)*
& HL1(L)*YLM(L,M)*GNT(L)
ENDIF
ENDDO
SUM_L=SUM_L*ATTKJ*4.D0*PI*IC
ELSE
SUM_L=ZEROC
ENDIF
C
IF(KLIN.EQ.JLIN) THEN
SM(KLIN,JLIN)=ONEC-TLK*
& SUM_L
ELSE
SM(KLIN,JLIN)=-TLK*SUM_L
ENDIF
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
LWORK=JLIN
C
C Inversion of the multiple scattering matrix MS and
C multiplication by T
C
CALL ZGETRF(JLIN,JLIN,SM,LINMAX*NATCLU_M,IPIV,INFO1)
IF(INFO1.NE.0) THEN
WRITE(6,*) ' ---> INFO1 =',INFO1
ELSE
CALL ZGETRI(JLIN,SM,LINMAX*NATCLU_M,IPIV,WORK,LWORK,INFO)
IF(INFO.NE.0) THEN
WRITE(6,*) ' ---> WORK(1),INFO =',WORK(1),INFO
ENDIF
ENDIF
C
C Storage of the Tau matrix
C
LIN=0
C
JLIN=0
DO JTYP=1,N_PROT
NBTYPJ=NATYP(JTYP)
LMJ=LMAX(JTYP,JE)
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
C
KLIN=0
DO KTYP=1,N_PROT
NBTYPK=NATYP(KTYP)
LMK=LMAX(KTYP,JE)
DO KNUM=1,NBTYPK
KATL=NCORR(KNUM,KTYP)
C
DO LJ=0,LMJ
ILJ=LJ*LJ+LJ+1
TLJ=DCMPLX(TL(LJ,1,JTYP,JE))
DO MJ=-LJ,LJ
INDJ=ILJ+MJ
JLIN=JLIN+1
C
DO LK=0,LMK
ILK=LK*LK+LK+1
DO MK=-LK,LK
INDK=ILK+MK
KLIN=KLIN+1
LIN=LIN+1
TAU(LIN)=CMPLX(SM(KLIN,JLIN)*TLJ)
ENDDO
ENDDO
KLIN=KLIN-INDK
C
ENDDO
ENDDO
KLIN=KLIN+INDK
JLIN=JLIN-INDJ
C
ENDDO
ENDDO
JLIN=JLIN+INDJ
C
ENDDO
ENDDO
C
RETURN
C
END

6809
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/lapack_inv.f Normal file
View File

File diff suppressed because it is too large Load Diff

915
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/leddif_mi.f Normal file
View File

@ -0,0 +1,915 @@
C
C
C=======================================================================
C
SUBROUTINE LEDDIF_MI(NPLAN,VAL,ZEM,IPHA,NAT2,COORD,NATYP,RHOK,
&NATCLU,NFICHLEC,JFICH,NP)
C
C This subroutine computes the LEED formula in the spin-independent
C case.
C
C The calculation is performed using a matrix inversion for the
C expression of the scattering path operator
C
C The matrix inversion is performed using the LAPACK inversion
C routines for a general complex matrix
C
C Last modified : 26 Apr 2013
C
C INCLUDE 'spec.inc'
USE DIM_MOD
USE ALGORITHM_MOD
USE AMPLI_MOD
USE APPROX_MOD
USE COOR_MOD, NTCLU => NATCLU, NTP => NATYP
USE DEBWAL_MOD
USE DIRECT_MOD, RTHETA => RTHEXT
USE EXTREM_MOD
USE FIXSCAN_MOD
USE INFILES_MOD
USE INUNITS_MOD
USE INIT_L_MOD
USE INIT_J_MOD
USE LIMAMA_MOD
USE MOYEN_MOD
USE OUTFILES_MOD
USE OUTUNITS_MOD
USE PARCAL_MOD
USE RESEAU_MOD
USE SPIN_MOD
USE TESTPB_MOD
USE TESTS_MOD
USE TRANS_MOD
USE TYPCAL_MOD
USE TYPEM_MOD
USE TYPEXP_MOD
USE VALIN_MOD
USE VALIN_AV_MOD
USE VALFIN_MOD
C
REAL BEAM(3),AXE(3),EPS(3),DIRBEAM(3),BEAMDIR(3),EMET(3)
C
COMPLEX IC,ONEC,ZEROC,COEF
COMPLEX TLT(0:NT_M,4,NATM,NE_M)
COMPLEX TAU(LINMAX*LINMAX*NATCLU_M*NATCLU_M)
COMPLEX YLMI(LINMAX)
COMPLEX YLMJ(LINMAX)
COMPLEX SLJDIF,SJDIF
COMPLEX SLIDIF,SLIDIR,SIDIF,SIDIR
COMPLEX RHOK(NE_M,NATM,0:18,2,NSPIN2_M),RD
COMPLEX ATT_MI,ATT_MI2,ATT_MJ
C
DIMENSION VAL(NATCLU_M),NATYP(NATM)
DIMENSION R_L(9),COORD(3,NATCLU_M)
C
C
C
CHARACTER*7 STAT
CHARACTER*13 OUTDATA1,OUTDATA2
C
C
CHARACTER*24 OUTFILE
CHARACTER*24 AMPFILE
C
C
DATA PI,PIS180,CONV /3.141593,0.017453,0.512314/
DATA FINSTRUC,CVECT,SMALL /0.007297,1.0,0.0001/
C
ALGO1='MI'
ALGO2=' '
ALGO3=' '
ALGO4=' '
C
I_DIR=0
NSET=1
JEL=1
JEMET=1
OUTDATA1='CROSS-SECTION'
IF(I_AMP.EQ.1) THEN
I_MI=1
OUTDATA2='MS AMPLITUDES'
ELSE
I_MI=0
ENDIF
C
C The first atom in the list taken as the origin
C
EMET(1)=SYM_AT(1,1)
EMET(2)=SYM_AT(2,1)
EMET(3)=SYM_AT(3,1)
C
IF(SPECTRO.EQ.'LED') THEN
IOUT=IUO2
OUTFILE=OUTFILE2
STAT='UNKNOWN'
IF(I_MI.EQ.1) THEN
IOUT2=IUSCR2+1
N_DOT=1
DO J_CHAR=1,24
IF(OUTFILE(J_CHAR:J_CHAR).EQ.'.') GOTO 888
N_DOT=N_DOT+1
ENDDO
888 CONTINUE
AMPFILE=OUTFILE(1:N_DOT)//'amp'
OPEN(UNIT=IOUT2, FILE=AMPFILE, STATUS=STAT)
ENDIF
ENDIF
C
C Position of the initial beam when the analyzer is along the z axis :
C (X_BEAM_Z,Y_BEAM_Z,Z_BEAM_Z)
C
RTHBEAM=THBEAM*PIS180
RPHBEAM=PHBEAM*PIS180
X_BEAM_Z=SIN(RTHBEAM)*COS(RPHBEAM)
Y_BEAM_Z=SIN(RTHBEAM)*SIN(RPHBEAM)
Z_BEAM_Z=COS(RTHBEAM)
C
IF(IMOD.EQ.0) THEN
C
C The analyzer is rotated
C
DIRBEAM(1)=X_BEAM_Z
DIRBEAM(2)=Y_BEAM_Z
DIRBEAM(3)=Z_BEAM_Z
ELSE
C
C The sample is rotated ---> beam and analyzer rotated
C
IF(I_EXT.EQ.0) THEN
RTH0=THETA0*PIS180
RPH0=PHI0*PIS180
RTH=RTH0
RPH=RPH0
C
C R_L is the rotation matrix from 0z to (THETA0,PHI0) expressed as
C a function of the Euler angles ALPHA=PHI0, BETA=THETA0, GAMMA=-PHI0
C It is stored as (1 2 3)
C (4 5 6)
C (7 8 9)
C
R_L(1)=COS(RTH0)*COS(RPH0)*COS(RPH0)+SIN(RPH0)*SIN(RPH0)
R_L(2)=COS(RTH0)*SIN(RPH0)*COS(RPH0)-SIN(RPH0)*COS(RPH0)
R_L(3)=SIN(RTH0)*COS(RPH0)
R_L(4)=COS(RTH0)*SIN(RPH0)*COS(RPH0)-SIN(RPH0)*COS(RPH0)
R_L(5)=COS(RTH0)*SIN(RPH0)*SIN(RPH0)+COS(RPH0)*COS(RPH0)
R_L(6)=SIN(RTH0)*SIN(RPH0)
R_L(7)=-SIN(RTH0)*COS(RPH0)
R_L(8)=-SIN(RTH0)*SIN(RPH0)
R_L(9)=COS(RTH0)
C
C Position of the beam when the analyzer is along (THETA0,PHI0) : BEAM(3)
C
BEAM(1)=X_BEAM_Z*R_L(1)+Y_BEAM_Z*R_L(2)+Z_BEAM_Z*R_L(3)
BEAM(2)=X_BEAM_Z*R_L(4)+Y_BEAM_Z*R_L(5)+Z_BEAM_Z*R_L(6)
BEAM(3)=X_BEAM_Z*R_L(7)+Y_BEAM_Z*R_L(8)+Z_BEAM_Z*R_L(9)
C
ENDIF
ENDIF
C
IC=(0.,1.)
ONEC=(1.,0.)
ZEROC=(0.,0.)
ATTSJ=1.
ATTSI=1.
ZSURF=VAL(1)
C
IF((ISOM.EQ.0).OR.(JFICH.EQ.1)) THEN
OPEN(UNIT=IOUT, FILE=OUTFILE, STATUS=STAT)
ENDIF
C
C Writing the headers in the output file
C
CALL HEADERS(IOUT)
C
IF((ISOM.EQ.0).OR.((ISOM.GT.0).AND.(JFICH.EQ.1))) THEN
WRITE(IOUT,12) SPECTRO,OUTDATA1
WRITE(IOUT,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,
& IE,IPH_1,I_EXT
IF(I_MI.EQ.1) THEN
WRITE(IOUT2,12) SPECTRO,OUTDATA2
WRITE(IOUT2,12) STEREO
WRITE(IOUT2,19) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,
& ITHETA,IE,IPH_1,I_EXT
WRITE(IOUT2,20) PHI0,THETA0,PHI1,THETA1,NONVOL(1)
ENDIF
ENDIF
C
IF(ISOM.EQ.0) THEN
WRITE(IOUT,79) NPLAN,NEMET,NTHETA,NPHI,NE
IF(I_MI.EQ.1) THEN
WRITE(IOUT2,79) NPLAN,NEMET,NTHETA,NPHI,NE
ENDIF
ELSEIF((ISOM.NE.0).AND.(JFICH.EQ.1)) THEN
WRITE(IOUT,11) NTHETA,NPHI,NE
IF(I_MI.EQ.1) THEN
WRITE(IOUT2,11) NTHETA,NPHI,NE
ENDIF
ENDIF
IJK=0
C
C Loop over the planes
C
DO JPLAN=1,NPLAN
Z=VAL(JPLAN)
IF((IPHA.EQ.1).OR.(IPHA.EQ.2)) THEN
DZZEM=ABS(Z-ZEM)
IF(DZZEM.LT.SMALL) GOTO 10
GOTO 1
ENDIF
10 CONTINUE
C
IF(ISOM.EQ.1) NP=JPLAN
C
C Loop over the energies
C
DO JE=1,NE
FMIN(0)=1.
FMAX(0)=1.
IF(NE.GT.1) THEN
ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1)
ELSEIF(NE.EQ.1) THEN
ECIN=E0
ENDIF
IF(I_TEST.NE.1) THEN
CFM=16.*PI*PI
ELSE
CFM=1.
ENDIF
CALL LPM(ECIN,XLPM,*6)
XLPM1=XLPM/A
IF(IPRINT.GT.0) WRITE(IUO1,56) A,XLPM1
IF((IPRINT.GT.0).AND.(IBAS.EQ.1)) THEN
IF(I_TEST.NE.2) WRITE(IUO1,57) COUPUR
ENDIF
IF(ITL.EQ.0) THEN
VK(JE)=SQRT(ECIN+VINT)*CONV*A*(1.,0.)
VK2(JE)=CABS(VK(JE)*VK(JE))
ENDIF
GAMMA=1./(2.*XLPM1)
IF(IPOTC.EQ.0) THEN
VK(JE)=VK(JE)+IC*GAMMA
ENDIF
IF(I_TEST.NE.1) THEN
VKR=REAL(VK(JE))
ELSE
VKR=1.
ENDIF
IF(I_MI.EQ.1) THEN
WRITE(IOUT2,21) ECIN,VKR*CFM
ENDIF
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) THEN
IF(IDCM.GE.1) WRITE(IUO1,22)
DO JAT=1,N_PROT
IF(IDCM.EQ.0) THEN
XK2UJ2=VK2(JE)*UJ2(JAT)
ELSE
XK2UJ2=VK2(JE)*UJ_SQ(JAT)
WRITE(IUO1,23) JAT,UJ_SQ(JAT)*A*A
ENDIF
CALL DWSPH(JAT,JE,XK2UJ2,TLT,ISPEED)
DO LAT=0,LMAX(JAT,JE)
TL(LAT,1,JAT,JE)=TLT(LAT,1,JAT,JE)
ENDDO
ENDDO
ENDIF
IF(ABS(I_EXT).GE.1) THEN
OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
READ(IUI6,13) I_DIR,NSET,N_DUM1
READ(IUI6,14) I_DUM1,N_DUM2,N_DUM3
ENDIF
C
C Largest angular momenum value (used to compute
C the spherical harmonics)
C
LM_MAX=0
DO JTYP=1,N_PROT
LMJ=LMAX(JTYP,JE)
LM_MAX=MAX(LM_MAX,LMJ)
ENDDO
C
C Initialization of TAU(LIN)
C
LIN=0
DO JTYP=1,N_PROT
NBTYPJ=NATYP(JTYP)
LMJ=LMAX(JTYP,JE)
DO JNUM=1,NBTYPJ
DO LJ=0,LMJ
DO MJ=-LJ,LJ
DO ITYP=1,N_PROT
NBTYPI=NATYP(ITYP)
LMI=LMAX(ITYP,JE)
DO INUM=1,NBTYPI
DO LI=0,LMI
DO MI=-LI,LI
LIN=LIN+1
TAU(LIN)=ZEROC
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C
C Matrix inversion for the calculation of TAU
C
IF(I_TEST.EQ.2) GOTO 666
CALL INV_MAT_MS(JE,TAU)
666 CONTINUE
C
C Calculation of the LEED formula
C
C
C Loop over the 'fixed' angle
C
15 DO J_FIXED=1,N_FIXED
IF(N_FIXED.GT.1) THEN
IF(I_EXT.EQ.0) THEN
FIX_STEP=(FIX1-FIX0)/FLOAT(N_FIXED-1)
XINCRF=FLOAT(J_FIXED-1)*FIX_STEP
ELSE
XINCRF=0.
ENDIF
ELSEIF(N_FIXED.EQ.1) THEN
XINCRF=0.
ENDIF
IF(ABS(I_EXT).GE.1) THEN
READ(IUI6,86) JSET,JLINE,THD,PHD
IF(I_EXT.EQ.-1) BACKSPACE IUI6
THETA0=THD
PHI0=PHD
ENDIF
IF(IPH_1.EQ.1) THEN
IF(I_EXT.EQ.0) THEN
DPHI=PHI0+XINCRF
ELSE
DPHI=PHD
ENDIF
RPHI=DPHI*PIS180
IF(IPRINT.GT.0) WRITE(IUO1,66) DPHI
ELSE
ISAUT=0
IF(I_EXT.EQ.0) THEN
DTHETA=THETA0+XINCRF
ELSE
DTHETA=THD
ENDIF
RTHETA=DTHETA*PIS180
IF(ABS(DTHETA).GT.90.) ISAUT=ISAUT+1
IF(I_EXT.GE.1) ISAUT=0
IF(I_TEST.EQ.2) ISAUT=0
IF(ISAUT.GT.0) GOTO 8
IF(IPRINT.GT.0) WRITE(IUO1,65) DTHETA
IF((IPRINT.GT.0).AND.(I_TEST.NE.2)) WRITE(IUO1,59)
IF((IPRINT.EQ.1).AND.(I_TEST.NE.2)) WRITE(IUO1,60)
C
C THETA-dependent number of PHI points for stereographic
C representation (to obtain a uniform sampling density).
C (Courtesy of J. Osterwalder - University of Zurich)
C
IF(STEREO.EQ.'YES') THEN
N_SCAN=INT((SCAN1-SCAN0)*SIN(RTHETA)/FIX_STEP+
& SMALL)+1
ENDIF
C
ENDIF
IF((N_FIXED.GT.1).AND.(IMOD.EQ.1)) THEN
C
C When there are several sets of scans (N_FIXED > 1),
C the initial position BEAM of the beam is recalculated
C for each initial position (RTH,RPH) of the analyzer
C
IF(IPH_1.EQ.1) THEN
RTH=THETA0*PIS180
RPH=RPHI
ELSE
RTH=RTHETA
RPH=PHI0*PIS180
ENDIF
C
R_L(1)=COS(RTH)*COS(RPH)
R_L(2)=-SIN(RPH)
R_L(3)=SIN(RTH)*COS(RPH)
R_L(4)=COS(RTH)*SIN(RPH)
R_L(5)=COS(RPH)
R_L(6)=SIN(RTH)*SIN(RPH)
R_L(7)=-SIN(RTH)
R_L(8)=0.
R_L(9)=COS(RTH)
C
BEAM(1)=X_BEAM_Z*R_L(1)+Y_BEAM_Z*R_L(2)+Z_BEAM_Z*R_L(3)
BEAM(2)=X_BEAM_Z*R_L(4)+Y_BEAM_Z*R_L(5)+Z_BEAM_Z*R_L(6)
BEAM(3)=X_BEAM_Z*R_L(7)+Y_BEAM_Z*R_L(8)+Z_BEAM_Z*R_L(9)
ENDIF
C
C Loop over the scanned angle
C
DO J_SCAN=1,N_SCAN
IF(N_SCAN.GT.1) THEN
XINCRS=FLOAT(J_SCAN-1)*(SCAN1-SCAN0)/FLOAT(N_SCAN-1)
ELSEIF(N_SCAN.EQ.1) THEN
XINCRS=0.
ENDIF
IF(I_EXT.EQ.-1) THEN
READ(IUI6,86) JSET,JLINE,THD,PHD
BACKSPACE IUI6
ENDIF
IF(IPH_1.EQ.1) THEN
ISAUT=0
IF(I_EXT.EQ.0) THEN
DTHETA=THETA0+XINCRS
ELSE
DTHETA=THD
ENDIF
RTHETA=DTHETA*PIS180
IF(ABS(DTHETA).GT.90.) ISAUT=ISAUT+1
IF(I_EXT.GE.1) ISAUT=0
IF(I_TEST.EQ.2) ISAUT=0
IF(ISAUT.GT.0) GOTO 8
IF(IPRINT.GT.0) WRITE(IUO1,65) DTHETA
IF((IPRINT.GT.0).AND.(I_TEST.NE.2)) WRITE(IUO1,59)
IF((IPRINT.EQ.1).AND.(I_TEST.NE.2)) WRITE(IUO1,60)
ELSE
IF(I_EXT.EQ.0) THEN
DPHI=PHI0+XINCRS
ELSE
DPHI=PHD
ENDIF
RPHI=DPHI*PIS180
IF(IPRINT.GT.0) WRITE(IUO1,66) DPHI
ENDIF
C
C Loop over the sets of directions to average over (for gaussian average)
C
C
SSETDIF=0.
SSETDIR=0.
C
SSET2DIF=0.
SSET2DIR=0.
C
IF(I_EXT.EQ.-1) THEN
JREF=INT(NSET)/2+1
ELSE
JREF=1
ENDIF
C
DO J_SET=1,NSET
IF(I_EXT.EQ.-1) THEN
READ(IUI6,86) JSET,JLINE,THD,PHD,W
DTHETA=THD
DPHI=PHD
RTHETA=DTHETA*PIS180
RPHI=DPHI*PIS180
C
C Here, there are several sets of scans (NSET > 1), so
C the initial position BEAM of the beam must be
C recalculated for each initial position of the analyzer
C
RTH=TH_0(J_SET)*PIS180
RPH=PH_0(J_SET)*PIS180
C
IF(IMOD.EQ.1) THEN
R_L(1)=COS(RTH)*COS(RPH)
R_L(2)=-SIN(RPH)
R_L(3)=SIN(RTH)*COS(RPH)
R_L(4)=COS(RTH)*SIN(RPH)
R_L(5)=COS(RPH)
R_L(6)=SIN(RTH)*SIN(RPH)
R_L(7)=-SIN(RTH)
R_L(8)=0.
R_L(9)=COS(RTH)
C
BEAM(1)=X_BEAM_Z*R_L(1)+Y_BEAM_Z*R_L(2)+
& Z_BEAM_Z*R_L(3)
BEAM(2)=X_BEAM_Z*R_L(4)+Y_BEAM_Z*R_L(5)+
& Z_BEAM_Z*R_L(6)
BEAM(3)=X_BEAM_Z*R_L(7)+Y_BEAM_Z*R_L(8)+
& Z_BEAM_Z*R_L(9)
C
ENDIF
ELSE
W=1.
ENDIF
C
IF(I_EXT.EQ.-1) PRINT 89
C
CALL DIRAN(VINT,ECIN,JEL)
C
IF(J_SET.EQ.JREF) THEN
DTHETAP=DTHETA
DPHIP=DPHI
ENDIF
C
IF(I_EXT.EQ.-1) THEN
WRITE(IUO1,88) DTHETA,DPHI
ENDIF
C
C .......... Case IMOD=1 only ..........
C
C Calculation of the position of the beam when the analyzer is at
C (THETA,PHI). DIRBEAM is the direction of the beam and its initial
C value (at (THETA0,PHI0)) is BEAM. AXE is the direction of the theta
C rotation axis and EPS is defined so that (AXE,DIRBEAM,EPS) is a
C direct orthonormal basis. The transform of a vector R by a rotation
C of OMEGA about AXE is then given by
C
C R' = R COS(OMEGA) + (AXE.R)(1-COS(OMEGA)) AXE + (AXE^R) SIN(OMEGA)
C
C Here, DIRANA is the internal direction of the analyzer and ANADIR
C its external position
C
C Note that when the initial position of the analyzer is (RTH,RPH)
C which coincides with (RTH0,RPH0) only for the first fixed angle
C
IF(IMOD.EQ.1) THEN
IF(ITHETA.EQ.1) THEN
AXE(1)=-SIN(RPH)
AXE(2)=COS(RPH)
AXE(3)=0.
RANGLE=RTHETA-RTH
ELSEIF(IPHI.EQ.1) THEN
AXE(1)=0.
AXE(2)=0.
AXE(3)=1.
RANGLE=RPHI-RPH
ENDIF
CALL PRVECT(AXE,BEAM,EPS,CVECT)
PRS=PRSCAL(AXE,BEAM)
IF(J_SCAN.EQ.1) THEN
DIRBEAM(1)=BEAM(1)
DIRBEAM(2)=BEAM(2)
DIRBEAM(3)=BEAM(3)
ELSE
DIRBEAM(1)=BEAM(1)*COS(RANGLE)+PRS*(1.-COS(
& RANGLE))*AXE(1)+SIN(RANGLE)*EPS(1)
DIRBEAM(2)=BEAM(2)*COS(RANGLE)+PRS*(1.-COS(
& RANGLE))*AXE(2)+SIN(RANGLE)*EPS(2)
DIRBEAM(3)=BEAM(3)*COS(RANGLE)+PRS*(1.-COS(
& RANGLE))*AXE(3)+SIN(RANGLE)*EPS(3)
ENDIF
ENDIF
IF(DIRBEAM(3).GT.1.) DIRBEAM(3)=1.
IF(DIRBEAM(3).LT.-1.) DIRBEAM(3)=-1.
THETABEAM=ACOS(DIRBEAM(3))
IF(I_TEST.EQ.2) THETABEAM=-THETABEAM
COEF=DIRBEAM(1)+IC*DIRBEAM(2)
CALL ARCSIN(COEF,DIRBEAM(3),PHIBEAM)
C
C Internal direction of the incoming beam BEAMDIR
C (DIRBEAM is the external direction)
C
CALL REFRAC(VINT,ECIN,THETABEAM,BEAMTHETA)
BEAMDIR(1)=SIN(BEAMTHETA)*COS(PHIBEAM)
BEAMDIR(2)=SIN(BEAMTHETA)*SIN(PHIBEAM)
BEAMDIR(3)=COS(BEAMTHETA)
C
CALL HARSPH3(NL_M,BEAMTHETA,-PHIBEAM,YLMI,LM_MAX)
C
ANABEAM=ANADIR(1,1)*DIRBEAM(1) + ANADIR(2,1)*DIRBEAM(2)
& +ANADIR(3,1)*DIRBEAM(3)
C
IF(IPRINT.GT.0) THEN
WRITE(IUO1,63) (DIRANA(J,1),J=1,3),(BEAMDIR(K),
& K=1,3),ANABEAM
ENDIF
IF(I_EXT.EQ.-1) PRINT 89
C
SRDIF=0.
SRDIR=0.
C
C Loop over the different directions of the analyzer contained in a cone
C
DO JDIR=1,NDIR
SIDIF=ZEROC
SIDIR=ZEROC
CALL HARSPH3(NL_M,THETAR(JDIR),PHIR(JDIR),YLMJ,
& LM_MAX)
C
C Loop over the first atom I encountered by the electron beam
C when entering the solid
C
LIN=0
DO ITYP=1,N_PROT
NBTYPI=NATYP(ITYP)
LMI=LMAX(ITYP,JE)
INDI_M=(LMI+1)*(LMI+1)
DO INUM=1,NBTYPI
IATL=NCORR(INUM,ITYP)
XOI=SYM_AT(1,IATL)-EMET(1)
YOI=SYM_AT(2,IATL)-EMET(2)
ZOI=SYM_AT(3,IATL)-EMET(3)
ROI=SQRT(XOI*XOI+YOI*YOI+ZOI*ZOI)
ZSURFI=VAL(1)-SYM_AT(3,IATL)
IF(IATTS.EQ.1) THEN
ATTSI=EXP(-ZSURFI*GAMMA/COS(BEAMTHETA))
ENDIF
IF(ROI.GT.SMALL) THEN
CSTHIR=(XOI*BEAMDIR(1)+YOI*BEAMDIR(2)+ZOI*
& BEAMDIR(3))/ROI
CTROIS1=ZOI/ROI
CSTHIR2=(XOI*(DIRANA(1,JDIR)-BEAMDIR(1))+YOI*
& (DIRANA(2,JDIR)-BEAMDIR(2))+ZOI*(DIRANA(3,JDIR)-
& BEAMDIR(3)))/ROI
ELSE
CSTHIR=0.
CTROIS1=0.
CSTHIR2=0.
ENDIF
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 78
IF(CTROIS1.GT.1.) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
IF(IDCM.GE.1) THEN
UJ2(ITYP)=UJ_SQ(ITYP)
ENDIF
IF(ABS(ZSURFI).LE.SMALL) THEN
IF(ABS(CSTHIR-1.).GT.SMALL) THEN
CSKZ2I=(BEAMDIR(3)-CTROIS1)*(BEAMDIR(3)-
& CTROIS1)/(2.-2.*CSTHIR)
ELSE
CSKZ2I=1.
ENDIF
UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.))
ELSE
UII=UJ2(ITYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UI2=VK2(JE)*UII
CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED)
ENDIF
78 IF(IDWSPH.EQ.1) THEN
DWTER=1.
DWTER2=1.
ELSE
DWTER=EXP(-VK2(JE)*UII*(1.-CSTHIR))
DWTER2=EXP(-VK2(JE)*UII*(1.-CSTHIR2))
ENDIF
ATT_MI=ATTSI*DWTER*CEXP(IC*VK(JE)*ROI*CSTHIR)
ATT_MI2=ATTSI*DWTER2*CEXP(-IC*VK(JE)*ROI*CSTHIR2)
C
C Kinematic term
C
SLIDIR=ZEROC
DO LI=0,LMI
ILI=LI*LI+LI+1
DO MI=-LI,LI
INDI=ILI+MI
SLIDIR=SLIDIR+TL(LI,1,ITYP,JE)*YLMJ(INDI)
& *YLMI(INDI)
ENDDO
ENDDO
C
C Loop over the last atom J encountered by the electron beam
C when exiting the solid
C
SJDIF=ZEROC
DO JTYP=1,N_PROT
NBTYPJ=NATYP(JTYP)
LMJ=LMAX(JTYP,JE)
INDJ_M=(LMJ+1)*(LMJ+1)
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
XOJ=SYM_AT(1,JATL)-EMET(1)
YOJ=SYM_AT(2,JATL)-EMET(2)
ZOJ=SYM_AT(3,JATL)-EMET(3)
ROJ=SQRT(XOJ*XOJ+YOJ*YOJ+ZOJ*ZOJ)
ZSURFJ=VAL(1)-SYM_AT(3,JATL)
IF(IATTS.EQ.1) THEN
ATTSJ=EXP(-ZSURFJ*GAMMA/DIRANA(3,JDIR))
ENDIF
IF(ROJ.GT.SMALL) THEN
CSTHJR=(XOJ*DIRANA(1,JDIR)+YOJ*DIRANA(2,
& JDIR)+ZOJ*DIRANA(3,JDIR))/ROJ
CTROIS1=ZOJ/ROJ
ELSE
CSTHJR=0.
CTROIS1=0.
ENDIF
IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 90
IF(CTROIS1.GT.1.) THEN
CTROIS1=1.
ELSEIF(CTROIS1.LT.-1.) THEN
CTROIS1=-1.
ENDIF
IF(IDCM.EQ.1) THEN
UJ2(JTYP)=UJ_SQ(JTYP)
ENDIF
IF(ABS(ZSURFJ).LE.SMALL) THEN
IF(ABS(CSTHJR-1.).GT.SMALL) THEN
CSKZ2J=(DIRANA(3,JDIR)-CTROIS1)*(
& DIRANA(3,JDIR)-CTROIS1)/(2.-2.*CSTHJR)
ELSE
CSKZ2J=1.
ENDIF
UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.))
ELSE
UJJ=UJ2(JTYP)
ENDIF
IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN
XK2UJ2=VK2(JE)*UJJ
CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED)
ENDIF
90 IF(IDWSPH.EQ.1) THEN
DWTER=1.
ELSE
DWTER=EXP(-VK2(JE)*UJJ*(1.-CSTHJR))
ENDIF
ATT_MJ=ATTSJ*DWTER*CEXP(-IC*VK(JE)*ROJ*
& CSTHJR)
C
C Loop over the angular momentum of atom I
C
SLIDIF=ZEROC
DO INDI=1,INDI_M
C
C Loop over the angular momentum of atom J
C
SLJDIF=ZEROC
DO INDJ=1,INDJ_M
LIN=LIN+1
SLJDIF=SLJDIF+YLMJ(INDJ)*TAU(LIN)
ENDDO
C
SLIDIF=SLIDIF+SLJDIF*YLMI(INDI)
ENDDO
C
C End of the loops over the last atom J
C
SJDIF=SJDIF+SLIDIF*ATT_MJ
C
ENDDO
ENDDO
SIDIF=SIDIF+SJDIF*ATT_MI
SIDIR=SIDIR+SLIDIR*ATT_MI2
C
C End of the loops over the first atom I
C
ENDDO
ENDDO
C
C Computing the square modulus
C
SRDIF=SRDIF+CABS(SIDIF)*CABS(SIDIF)
SRDIR=SRDIR+CABS(SIDIR)*CABS(SIDIR)
C
C End of the loop on the directions of the analyzer
C
ENDDO
C
SSETDIF=SSETDIF+SRDIF*CFM*W/NDIR
SSETDIR=SSETDIR+SRDIR*CFM*W/NDIR
IF(ICHKDIR.EQ.2) THEN
IF(JSET.EQ.JREF) THEN
SSET2DIF=SRDIF*CFM/NDIR
SSET2DIR=SRDIR*CFM/NDIR
ENDIF
ENDIF
C
C End of the loop on the set averaging
C
ENDDO
C
IF(ISOM.EQ.2) THEN
WRITE(IOUT,67) JPLAN,JFICH,DTHETAP,DPHIP,ECIN,
& SSETDIR,SSETDIF
IF(ICHKDIR.EQ.2) THEN
WRITE(IOUT,67) JPLAN,JFICH,DTHETAP,DPHIP,ECIN,
& SSET2DIR,SSET2DIF
ENDIF
ELSE
WRITE(IOUT,67) JPLAN,JEMET,DTHETAP,DPHIP,ECIN,
& SSETDIR,SSETDIF
IF(ICHKDIR.EQ.2) THEN
WRITE(IOUT,67) JPLAN,JEMET,DTHETAP,DPHIP,ECIN,
& SSET2DIR,SSET2DIF
ENDIF
ENDIF
C
C End of the loop on the scanned angle
C
ENDDO
C
8 CONTINUE
C
C End of the loop on the fixed angle
C
ENDDO
C
C End of the loop on the energy
C
CLOSE(IUI6)
ENDDO
C
3 CONTINUE
C
GO TO 1
5 IPLAN=JPLAN-1
IJK=IJK+1
IF((IJK.EQ.1).AND.(IPRINT.GT.0)) THEN
IF(I_TEST.NE.2) WRITE(IUO1,54) IPLAN
ENDIF
1 CONTINUE
C
C End of the loop on the planes
C
ENDDO
C
IF(ABS(I_EXT).GE.1) CLOSE(IUI6)
IF((ISOM.EQ.0).OR.(JFICH.EQ.NFICHLEC)) WRITE(IOUT,*)
IF(SPECTRO.EQ.'APC') CLOSE(IOUT)
IF(SPECTRO.EQ.'APC') GOTO 7
c IF(((NEMET.GT.1).OR.(NPLAN.GT.1)).AND.(ISOM.EQ.0)) THEN
IF(((NEMET.GT.1).OR.(NPLAN.GT.0)).AND.(ISOM.EQ.0)) THEN
NP=0
CALL TREAT_PHD(ISOM,NFICHLEC,JFICH,NP)
ENDIF
IF(I_EXT.EQ.2) THEN
CALL WEIGHT_SUM(ISOM,I_EXT,0,1)
ENDIF
GOTO 7
6 WRITE(IUO1,55)
C
9 FORMAT(9(2X,I1),2X,I2)
11 FORMAT(I4,2X,I4,2X,I4)
12 FORMAT(2X,A3,11X,A13)
13 FORMAT(6X,I1,1X,I3,2X,I4)
14 FORMAT(6X,I1,1X,I3,3X,I3)
19 FORMAT(2(2X,I1),1X,I2,6(2X,I1),2X,I2)
20 FORMAT(2(5X,F6.2,2X,F6.2),2X,I1)
21 FORMAT(10X,E12.6,3X,E12.6)
22 FORMAT(16X,'INTERNAL CALCULATION OF MEAN SQUARE DISPLACEMENTS',/,
&25X,' BY DEBYE UNCORRELATED MODEL:',/)
23 FORMAT(21X,'ATOM TYPE ',I5,' MSD = ',F8.6,' ANG**2')
51 FORMAT(/////,2X,'******* PLANE NUMBER ',I3,' DOES NOT CONTAIN ',
*'ANY ABSORBER OF TYPE ',I2,' *******')
52 FORMAT(/////,2X,'******* PLANE NUMBER ',I3,' POSITION OF ','THE
&ABSORBER : (',F6.3,',',F6.3,',',F6.3,') *******',/,2X,'******* ',
&19X,'THIS ABSORBER IS OF TYPE ',I2,20X,' *******')
53 FORMAT(//,2X,'ORDER ',I2,' TOTAL NUMBER OF PATHS : ',F15.1,/
&,10X,' EFFECTIVE NUMBER OF PATHS : ',F15.1,/,10X,' MINIMAL
&INTENSITY : ',E12.6,2X,'No OF THE PATH : ',F15.1,
& /,10X,' MAXIMAL INTENSITY : ',E12.6,2X,
&'No OF THE PATH : ',F15.1)
54 FORMAT(//,7X,'DUE TO THE SIZE OF THE CLUSTER, THE SUMMATION',
*' HAS BEEN TRUNCATED TO THE ',I2,' TH PLANE')
55 FORMAT(///,12X,' <<<<<<<<<< THIS VALUE OF ILPM IS NOT',
*'AVAILABLE >>>>>>>>>>')
56 FORMAT(4X,'LATTICE PARAMETER A = ',F6.3,' ANGSTROEMS',4X,
*'MEAN FREE PATH = ',F6.3,' * A',//)
57 FORMAT(25X,'CLUSTER RADIUS = ',F6.3,' *A')
58 FORMAT(//,2X,'ORDER ',I2,' TOTAL NUMBER OF PATHS : ',I10,/,
&10X,' EFFECTIVE NUMBER OF PATHS : ',I10, /,10X,'
& MINIMAL INTENSITY : ',E12.6,2X,'No OF THE PATH : ',I10,
& /,10X,' MAXIMAL INTENSITY : ',
&E12.6, 2X,'No OF THE PATH : ',I10)
59 FORMAT(//,15X,'THE SCATTERING DIRECTION IS GIVEN INSIDE ',
*'THE CRYSTAL')
60 FORMAT(7X,'THE POSITIONS OF THE ATOMS ARE GIVEN WITH RESPECT ',
*'TO THE ABSORBER')
63 FORMAT(///,4X,'.......... DIRECTION OF THE DETECTOR : (',
&F6.3,',',F6.3,',',F6.3, ') ..........',/,16X,'DIRECTION OF
&THE BEAM ', ' : (',F6.3,',',F6.3,',',F6.3,')',/,16X,
&'ANALYZER.BEAM : ',F7.4)
65 FORMAT(////,3X,'++++++++++++++++++',9X,
*'THETA = ',F6.2,' DEGREES',9X,'++++++++',
*'++++++++++',///)
66 FORMAT(////,3X,'++++++++++++++++++',9X,
*'PHI = ',F6.2,' DEGREES',9X,'++++++++++',
*'++++++++++',///)
67 FORMAT(2X,I3,2X,I2,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
68 FORMAT(10X,' CUT-OFF INTENSITY : ',E12.6)
69 FORMAT(9X,I2,2X,E12.6,7X,E12.6,1X,F6.3,1X,10(I3,2X))
70 FORMAT(2X,I2,2X,I10,7X,E12.6,2X,F6.3,7X,I2,7X,10(I3,2X))
71 FORMAT(//,1X,'JDIF',4X,'No OF THE PATH',2X,'INTENSITY',3X,
&'LENGTH',4X,'ABSORBER',2X,'ORDER OF THE SCATTERERS',/)
72 FORMAT(2X,I3,2X,I2,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,
&E12.6,2X,E12.6)
74 FORMAT(10X,'<===== NUMBER OF PATHS TOO LARGE FOR PRINTING ','====
&=>')
76 FORMAT(2X,I2,2X,E12.6,7X,E12.6,2X,F6.3,7X,I2,7X,10(I3,2X))
77 FORMAT(' ')
79 FORMAT(2X,I3,2X,I2,2X,I4,2X,I4,2X,I4)
80 FORMAT(///)
81 FORMAT(//,1X,'RANK',1X,'ORDER',4X,'No PATH',3X,'INTENSITY',3X,
&'LENGTH',4X,'ABS',3X,'ORDER OF THE SCATTERERS',/)
82 FORMAT(I3,4X,I2,1X,E12.6,3X,E12.6,2X,F6.3,4X,I2,4X,10(I3,1X))
83 FORMAT(I3,4X,I2,1X,I10,3X,E12.6,2X,F6.3,4X,I2,4X,10(I3,1X))
84 FORMAT(/////,18X,'THE ',I3,' MORE INTENSE PATHS BY DECREASING','
&ORDER :',/,24X,'(THE LENGTH IS GIVEN IN UNITS ','OF A)')
85 FORMAT(/////,25X,' PATHS USED IN THE CALCULATION : ',/,24X,'(THE
&LENGTH IS GIVEN IN UNITS OF A)')
86 FORMAT(2X,I3,1X,I4,5X,F8.3,3X,F8.3,3X,E12.6)
87 FORMAT(2X,I2,2X,I3,2X,I2,2X,I3,2X,I3,2X,I3,2X,I1,2X,I2,2X,I2,2X,
&E12.6,2X,E12.6,2X,E12.6,2X,E12.6)
88 FORMAT(/,19X,'TILTED THETA =',F6.2,5X,'TILTED PHI =', F6.2)
89 FORMAT(/,4X,'..........................................','.......
&..............................')
C
7 RETURN
C
END

21
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/main.f Normal file
View File

@ -0,0 +1,21 @@
SUBROUTINE RUN(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_,
& NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_,
& NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_,
& N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
USE DIM_MOD
IMPLICIT INTEGER (A-Z)
CF2PY INTEGER, INTENT(IN,COPY) :: NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_
CF2PY INTEGER, INTENT(IN,COPY) :: NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_
CF2PY INTEGER, INTENT(IN,COPY) :: NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_
CF2PY INTEGER, INTENT(IN,COPY) :: N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_
CALL ALLOCATION(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_,
& NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_,
& NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_,
& N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
CALL MAIN_LED_NS_MI()
CALL CLOSE_ALL_FILES()
END SUBROUTINE RUN

1648
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/main_led_ns_mi.f Normal file
View File

File diff suppressed because it is too large Load Diff

106
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/plotfd.f Normal file
View File

@ -0,0 +1,106 @@
C
C=======================================================================
C
SUBROUTINE PLOTFD(A,LMX,ITL,NL,NAT,NE)
C
C This routine prepares the output for a plot of the scattering factor
C
USE DIM_MOD
C
USE APPROX_MOD
USE FDIF_MOD
USE INIT_L_MOD , L => LI, I2 => INITL, I3 => NNL, I4 => LF1, I5 =>
& LF2, I10 => ISTEP_LF
USE INIT_J_MOD
USE OUTFILES_MOD
USE OUTUNITS_MOD
USE PARCAL_MOD , N3 => NPHI, N4 => NE, N5 => NTHETA, N6 => NEPS
USE TYPCAL_MOD , I7 => IFTHET, I8 => IMOD, I9 => IPOL, I12 => I_CP
&, I13 => I_EXT, I14 => I_TEST
USE VALIN_MOD , U1 => THLUM, U2 => PHILUM, U3 => ELUM, N7 => NONVO
&L
USE VALFIN_MOD
C
C
C
DIMENSION LMX(NATM,NE_M)
C
COMPLEX FSPH,VKE
C
C
C
DATA PI,CONV/3.141593,0.512314/
C
OPEN(UNIT=IUO3, FILE=OUTFILE3, STATUS='UNKNOWN')
IF(ISPHER.EQ.0) THEN
L=0
LMAX=0
ELSE
LMAX=L
ENDIF
PHITOT=360.
THTOT=360.*ITHETA*(1-IPHI)+180.*ITHETA*IPHI
NPHI=(NFTHET+1)*IPHI+(1-IPHI)
NTHT=(NFTHET+1)*ITHETA*(1-IPHI)+(NFTHET/2+1)*ITHETA*IPHI+
* (1-ITHETA)
NE=NFTHET*IE + (1-IE)
WRITE(IUO3,1) ISPHER,NL,NAT,L,NTHT,NPHI,NE,E0,EFIN
DO 10 JT=1,NTHT
DTHETA=THETA1+FLOAT(JT-1)*THTOT/FLOAT(MAX0(NTHT-1,1))
RTHETA=DTHETA*PI/180.
TEST=SIN(RTHETA)
IF(TEST.GE.0.) THEN
POZ=PI
EPS=1.
ELSE
POZ=0.
EPS=-1.
ENDIF
BETA=RTHETA*EPS
IF(ABS(TEST).LT.0.0001) THEN
NPHIM=1
ELSE
NPHIM=NPHI
ENDIF
DO 20 JP=1,NPHIM
DPHI=PHI1+FLOAT(JP-1)*PHITOT/FLOAT(MAX0(NPHI-1,1))
RPHI=DPHI*PI/180.
GAMMA=POZ-RPHI
DO 30 JE=1,NE
IF(NE.EQ.1) THEN
ECIN=E0
ELSE
ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1)
ENDIF
IF(ITL.EQ.0) VKE=SQRT(ECIN-ABS(VINT))*CONV*A*(1.,0.)
DO 40 JAT=1,NAT
IF(L.GT.LMX(JAT,JE)) GOTO 90
DO 50 M=-LMAX,LMAX
CALL FACDIF1(VKE,R1,R2,THETA0,PHI0,BETA,GAMMA,L,M,FSPH,J
&AT,JE,*60)
GOTO 70
60 WRITE(IUO1,80)
STOP
70 REFTH=REAL(FSPH)
XIMFTH=AIMAG(FSPH)
WRITE(IUO3,5) JE,JAT,L,M,REFTH,XIMFTH,DTHETA,DPHI,ECIN
50 CONTINUE
GOTO 40
90 WRITE(IUO1,100) JAT
STOP
40 CONTINUE
30 CONTINUE
20 CONTINUE
10 CONTINUE
CLOSE(IUO3)
1 FORMAT(5X,I1,2X,I2,2X,I4,2X,I2,2X,I3,2X,I3,2X,I3,2X,F8.2,2X,F8.2)
5 FORMAT(1X,I3,1X,I4,1X,I2,1X,I3,1X,F6.3,1X,F6.3,1X,F6.2,1X,F6.2,1X,
&F8.2)
80 FORMAT(15X,'<<<<< WRONG VALUE OF THETA0 : THE DENOMINATOR ','IS Z
&ERO >>>>>')
100 FORMAT(15X,'<<<<< THE VALUE OF L EST IS TOO LARGE FOR ATOM',' : '
&,I2,' >>>>>')
C
RETURN
C
END

995
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/process.py Normal file
View File

@ -0,0 +1,995 @@
# vim: set et ts=4 sw=4 fdm=indent:
# coding: utf-8
import re
import sys
import os
import textwrap
class Patterns(object):
col = '(?P<comment> |C|\*)'
col += '(?P<label>(?: |\d){1,5})'
col += '(?P<cont> |\d|&)'
typ = '(?P<type>'
typ += 'BYTE|'
typ += 'CHARACTER(?:\*\(\*\)|\*\d+)?|'
typ += 'COMPLEX(?:\*8|\*16|\*32)|'
typ += 'DOUBLE\s+(?:COMPLEX|PRECISION)|'
typ += 'INTEGER(?:\*2|\*4|\*8)?|'
typ += 'LOGICAL(?:\*1|\*2|\*4|\*8)?|'
typ += 'REAL(?:\*4|\*8|\*16)?|'
typ += 'AUTOMATIC|STATIC)'
nam = '[a-z][a-z0-9_]*'
dim = '(?:,?.*?(?::.*?)?)+'
axs = '[^:]+(:[^:]+)?'
class BaseInfo(object):
def __init__(self, **kwargs):
self._attrs = kwargs.keys()
for kw,val in kwargs.items():
setattr(self, '_' + kw, val)
@property
def info(self):
s = '=== {}:\n'.format(self.__class__.__name__)
for attr in self._attrs:
s += ' {}: {}\n'.format(attr, repr(getattr(self, attr)))
return s
def __repr__(self):
return '<{}>'.format(self.__class__.__name__)
class DimensionInfo(BaseInfo):
def __init__(self, **kwargs):
opts = {'rank': None,
'extents': None,
'variable': None}
opts.update(**kwargs)
BaseInfo.__init__(self, **opts)
@property
def rank(self):
return len(self.extents)
@property
def extents(self):
return self._extents
@property
def variable(self):
return self._variable
@variable.setter
def variable(self, value):
assert isinstance(value, VariableInfo)
self._variable = value
def __str__(self):
s = ''
for d in self.extents:
s += d[0]
if d[1] is not None:
s += ':' + d[1]
s += ','
s = s.strip(',')
return s
class VariableInfo(BaseInfo):
def __init__(self, **kwargs):
opts = {'name': None, 'type': None,
'dimension': None, 'subprogram': None}
opts.update(**kwargs)
BaseInfo.__init__(self, **opts)
@property
def name(self):
return self._name
@property
def type(self):
return self._type
@type.setter
def type(self, value):
self._type = value
@property
def dimension(self):
return self._dimension
@dimension.setter
def dimension(self, value):
self._dimension = value
@property
def subprogram(self):
return self._subprogram
@subprogram.setter
def subprogram(self, value):
self._subprogram = value
def __str__(self):
s = self.name
if self.dimension is not None:
s += '(' + str(self.dimension) + ')'
return s
class FileInfo(BaseInfo):
def __init__(self, **kwargs):
opts = {'filename': None, 'content': None, 'subprograms': None}
opts.update(**kwargs)
BaseInfo.__init__(self, **opts)
self._subprograms = find_subprograms('\n'.join(self.content))
for sp in self._subprograms:
sp.file = self
#self._subprograms = None
@property
def filename(self):
return os.path.abspath(self._filename)
@property
def content(self):
with open(self.filename, 'r') as fd:
lines = fd.readlines()
pat = re.compile(' (?:\d|&)\s*(.*)$')
c = []
for line in lines:
line = line.strip('\n')
m = pat.match(line)
if m:
c[-1] += m.group(1)
else:
c.append(line)
return c
@property
def subprograms(self):
return self._subprograms
def __str__(self):
return content2str(self.content)
class SubProgramInfo(BaseInfo):
def __init__(self, **kwargs):
opts = {'name': None, 'content': None, 'type': None, 'file': None,
'l0': None, 'l1': None, 'commons': None}
opts.update(**kwargs)
BaseInfo.__init__(self, **opts)
@property
def name(self):
return self._name
@property
def content(self):
return self._file.content[self.l0:self.l1]
@property
def type(self):
return self._type
@property
def file(self):
return self._file
@property
def l0(self):
return self._l0
@property
def l1(self):
return self._l1
@file.setter
def file(self, value):
self._file = value
@property
def commons(self):
self._commons = find_commons('\n'.join(self.content))
for c in self._commons:
c.subprogram = self
return self._commons
def __str__(self):
return content2str(self.content)
def __repr__(self):
s = "{}<{}>".format(self.name, self.type)
return s
class CommonBlockInfo(BaseInfo):
def __init__(self, **kwargs):
opts = {'name': None, 'content': None, 'subprogram': None, 'variables': None}
opts.update(**kwargs)
BaseInfo.__init__(self, **opts)
def __str__(self):
return content2str(self.content)
def __repr__(self):
s = "{}<{:d} variables>".format(self.name, len(self.variables))
return s
@property
def subprogram(self):
return self._subprogram
@subprogram.setter
def subprogram(self, value):
self._subprogram = value
@property
def name(self):
return self._name
@property
def content(self):
return self._content
@property
def variables(self):
# string to analyse
s = self.content[0]
m = re.match("^.*COMMON\s*/{}/\s*(.*)$".format(self.name),s, re.I)
self._variables = find_variables(m.groups()[0])
for v in self._variables:
v.subprogram = self.subprogram
# If the dimension of a variable is None, try to search if a DIMENSION
# statement exists in the subprogram
dim_defs = [] # list of variables defined in a DIMENSION statement
for line in self.subprogram.content:
m = re.match("^[^C]\s+DIMENSION\s+(.*)$", line, re.I)
if m is not None:
s = m.groups()[0]
var_list = find_variables(s)
for v in var_list:
dim_defs.append(v)
dim_defs = Variables(dim_defs)
# Now for each variable of the common, if there is no dimension, try to find it in
# the dim_defs list
for v in self._variables:
if v.dimension is None:
dim = dim_defs[v.name]
if dim is not None:
#print(self.subprogram.name)
#print(v.name, dim)
v.dimension = dim.dimension
#exit()
# if the type of the variable is None, try to find it in a declaration statement
type_defs = [] # list of variables defined with their type
for line in self.subprogram.content:
var_list = find_type(line)
if var_list is not None:
for v in var_list:
type_defs.append(v)
type_defs = Variables(type_defs)
# Now for each variable with no type, try to find it in the type_defs list
for v in self._variables:
if v.type is None:
typ = type_defs[v.name]
if typ is not None:
#print(self.subprogram.name)
#print(v.name, typ)
v.type = typ.type
#exit()
else:
if re.match("^[A-HO-Z].*", v.name, re.I):
v.type = "REAL"
else:
v.type = "INTEGER"
return self._variables
class _CommonBlockInfo(BaseInfo):
def __init__(self, **kwargs):
self.name = kwargs.get('name', None)
self.subprogram = kwargs.get('subprogram', None)
def __init__(self, **kwargs):
opts = {'name': None, 'content': None, 'type': None, 'file': None,
'l0': None, 'l1': None}
opts.update(**kwargs)
BaseInfo.__init__(self, **opts)
def __str__(self):
s = "Common block name: {}\n".format(self.name)
for var in self.variables:
s += str(var)
return s
def __repr__(self):
s = '{}'.format(self.name)
return s
def _find_variables(self):
content = self.content[0].rstrip()
pat = re.compile('^\s*COMMON\s+/{}/\s+(.*)$'.format(self.name), re.IGNORECASE)
m = pat.match(content)
var_loc = m.group(1)
p0 = re.compile(r'\(.*[^\)]$')
p1 = re.compile('^[a-zA-Z0-9_\*]*\)$')
var_list = []
var_loc_list = var_loc.split(',')
for i, _ in enumerate(var_loc_list):
_ = _.strip()
if i > 0:
if p0.search(var_list[-1]):
var_list[-1] += ',' + _
else:
var_list.append(_)
else:
var_list.append(_)
variables = []
for var in var_list:
m = re.match("([a-zA-Z0-9_]+)(\((.*)\))?", var)
var_name = m.group(1)
v = VariableInfo(name=var_name, subprogram=self.subprogram)
#dim_loc = m.groups()[-1]
#if dim_loc is not None:
# dim_list = dim_loc.split(',')
# v.dimension = dim_list
variables.append(v)
return variables
@property
def content(self):
pat = re.compile('\s+.*COMMON\s+/{}/.*$'.format(self.name), re.IGNORECASE)
for line in self.subprogram.content:
if pat.match(line):
return [line.rstrip(),]
@property
def variables(self):
variables_list = self._find_variables()
return Variables(variables_list)
class InfoList(object):
def __init__(self, elements):
self._elements = elements
for element in self._elements:
setattr(self, element.name, element)
def __getitem__(self, item):
if isinstance(item, str):
for element in self._elements:
if element.name == item:
return element
return None
elif isinstance(item, int):
return self._elements[item]
else:
raise NameError('Unable to retrieve item!')
def __str__(self):
s = "nb of {}: {:d}\n".format(self.__class__.__name__, len(self._elements))
#s += str([_.name for _ in self._elements])
for _ in self._elements:
s += _.info
return s
def __len__(self):
return len(self._elements)
class Subprograms(InfoList):
def __init__(self, *args, **kwargs):
InfoList.__init__(self, *args, **kwargs)
class Commons(InfoList):
def __init__(self, *args, **kwargs):
InfoList.__init__(self, *args, **kwargs)
class Variables(InfoList):
def __init__(self, *args, **kwargs):
InfoList.__init__(self, *args, **kwargs)
"""
class Subprograms(BaseInfo):
def __init__(self, subprograms_list):
self._subprograms = subprograms_list
for sp in self._subprograms:
setattr(self, sp.name, sp)
def __getitem__(self, item):
return self._subprograms[item]
def __str__(self):
s = "nb of subprograms: {:d}\n".format(len(self._subprograms))
return s
class Commons(BaseInfo):
def __init__(self, commons_list):
self._commons = commons_list
for cmn in self._commons:
setattr(self, cmn.name, cmn)
def __getitem__(self, item):
return self._commons[item]
def __str__(self):
s = "nb of commons: {:d}\n".format(len(self._commons))
return s
class Variables(BaseInfo):
"""
class VariableInfo2(BaseInfo):
def __init__(self, **kwargs):
self.name = kwargs.get('name', None)
#self.type = kwargs.get('type', None)
self.dimension = kwargs.get('dimension', None)
self.i = []
self.j = []
self.subprogram = kwargs.get('subprogram', None)
#self._find_type()
#self._find_dimension()
def __str__(self):
s = "Variable name: {}\n".format(self.name)
s += " type: {}\n".format(self.type)
s += " dimension: {}\n".format(self.dimension)
return s
def _find_implicit(self):
content = self.subprogram.content
pat = re.compile('^\s*IMPLICIT\s+', re.IGNORECASE)
def _find_type(self):
content = self.subprogram.content
pat = re.compile('^\s+((?:INTEGER|REAL|DOUBLE PRECISION|COMPLEX|LOGICAL|CHARACTER)\S*).*{}[\(,]?.*$'.format(self.name), re.IGNORECASE)
for line in content:
m = pat.match(line)
print(line)
if m:
return m.group(1).strip()
return None
def _find_dimension(self):
content = self.subprogram.content
dimension = None
#pat = re.compile('^\s+DIMENSION.*{}\(([^\(])\).*$'.format(self.name),re.IGNORECASE)
pat = re.compile('^\s+DIMENSION.*{}\((.*?)\).*$'.format(self.name),re.IGNORECASE)
for line in content:
m = pat.match(line)
if m:
print(line)
#dimension = m.group(1).strip()
dimension = m.group(1)
print(dimension)
if dimension is not None:
print('Variable: {}, dimension: {}'.format(self.name, dimension))
@property
def type(self):
t = self._find_type()
return t
###############################################################################
# UTILITY FUNCTIONS
###############################################################################
def splitline_(line, width=72):
result = []
i = 0
j = len(line)
ll = line[i:j]
if len(ll) > width:
s = ''
for dec in range(8):
s += '{:d} '.format(dec)
print(s)
print('0123456789'*7+'012')
#print(line)
while len(ll) > width:
breaks = [_.end() for _ in re.finditer('[ ,]', ll)]
print(ll)
print(breaks)
for ij, j in enumerate(breaks):
tmp = ll[i:j]
#print(breaks,i,j,ij, tmp, len(tmp))
if len(tmp) <= width and ij < len(breaks)-1:
continue
else:
_ = ll[i:breaks[ij-1]]
result.append(_)
i = len(_)
if i <= 6:
print(j, _)
raise NameError('Impossible to cut line at breaks')
ll = ' &' + ll[i:]
i = 0
break
result.append(ll)
print(ll)
return result
def splitline(line, width=72):
if len(line) == 0:
return []
if line[0].upper() == 'C':
return [line,]
head = line[:6]
L = line[6:] # the working line
# find the indentation
m = re.search('^\s*', L)
indent = L[m.start():m.end()]
# and define the true width to work with
W = width - 6 - len(indent)
def rule():
s = ''
for dec in range(20):
s += '{:d} '.format(dec)
print(s)
print('0123456789'*20)
# find all places to break the line
breaks = [_.end() for _ in re.finditer('[ ,()+-/*=:]', L)]
# split at breaks
indices = [0,] + breaks + [len(L),]
indices = zip(indices[:-1], indices[1:])
splitted_line = [L[a:b] for a,b in indices]
# iterate over each element and add it to the previous one if length is < max
chain = [splitted_line[0],]
for element in splitted_line[1:]:
l1 = len(chain[-1])
l2 = len(element)
if l1+l2 < W:
chain[-1] = chain[-1] + element
else:
chain.append(element)
# restore the head of the line
chain[0] = head + chain[0]
# add the & symbol
for i in range(1,len(chain)):
chain[i] = " &" + indent + chain[i]
# final check
for element in chain:
if len(element) > width:
rule()
print(f"{line}")
print(f"breaks at = {breaks}")
print(chain)
rule()
print(element)
raise NameError(f"Unable to split line!")
return chain
def content2str(content):
new_content = []
for index, line in enumerate(content):
#print(f'{index:>5d}#{line}')
multilines = splitline(line.rstrip(), width=72)
#print(multilines)
new_content += multilines
return '\n'.join(new_content)
def split_at_comma(string):
"""
"""
# remove all spaces from the string
line0 = string.replace(' ', '')
line = line0
# define some patterns
pat0 = re.compile('\(([^\(\)]*)\)', re.I)
# remove nested blocks in ()'s and replace them with
# hash signs (#) to make the treatment easier
while True:
M = list(pat0.finditer(line))
if len(M) == 0: break
for m in M:
i,j = m.start(), m.end()
line = line[:i] + '#'*(j-i) + line[j:]
# now get indices of ','
indices = [_.start() for _ in re.finditer(',', line)]
indices = zip([-1,] + indices, indices + [len(line),])
# create the list of
elements = [line0[i+1:j] for i,j in indices]
return elements
def find_dimension(string):
"""
Finds the components of a dimension declaration.
:param string: The argument of a dimension declaration
:type string: str
:return: A DimensionInfo object.
:rtype: DimensionInfo
:Example:
>>> dim = find_dimension('I,J,-3:2')
>>> print(dim)
>>> (3, ())
"""
# define some patterns
pat0 = re.compile('([^:]+):?([^:]+)?', re.I)
# create the list of axes
axl = split_at_comma(string)
# get the extents
extents = []
for ax in axl:
m = pat0.match(ax)
extents.append(m.groups())
return DimensionInfo(extents=extents)
def find_variables(string):
"""
Finds the name and dimension of variables in a comma separated
list of variables.
:param string: The comma separated variables declaration
:type string: str
:return: A Variables object.
:rtype: Variables
:Example:
>>> variables = find_variables('ONE, TWO(3,3)')
>>> print(variables)
>>> nb of Variables: 2
>>> === VariableInfo:
>>> name: 'ONE'
>>> type: None
>>> dimension: None
>>> subprogram: None
>>> === VariableInfo:
>>> name: 'TWO'
>>> type: None
>>> dimension: <DimensionInfo>
>>> subprogram: None
"""
# create the list of variables
var_list = []
variables = split_at_comma(string)
pat0 = re.compile('({})(?:\((.*)\))?'.format(Patterns.nam), re.I)
for var in variables:
# extract the variable's name and dimension if any
m = pat0.match(var)
name = m.group(1)
if m.group(2) is not None:
dimension = find_dimension(m.group(2))
else:
dimension = None
variable = VariableInfo(name=name, dimension=dimension)
if isinstance(dimension, DimensionInfo):
dimension.variable = variable
var_list.append(variable)
return Variables(var_list)
def find_subprograms(string):
lines = string.split('\n')
subprograms = []
for iline, line in enumerate(lines):
patterns = [('SUBROUTINE', re.compile("\s*SUBROUTINE\s+(\w+)\(?.*")),
('FUNCTION', re.compile("\s*.*FUNCTION\s+(\w+)\(?.*")),
('PROGRAM', re.compile("\s*PROGRAM\s+(\w+).*"))]
for t, pat in patterns:
m = pat.match(line)
if m is not None:
subprog = SubProgramInfo(type=t,
name=m.group(1),
l0=iline)
subprograms.append(subprog)
for i, subprog in enumerate(subprograms):
if i < len(subprograms) - 1:
subprog._l1 = subprograms[i+1].l0
else:
subprog._l1 = -1
return Subprograms(subprograms)
def find_commons(string):
pat = re.compile("^\s+COMMON\s*/([a-zA-Z0-9_]+)/(.*)$")
commons = []
for line in string.split('\n'):
# extract the name of the common block
m = pat.match(line)
if m is not None:
# name
name = m.group(1)
c = CommonBlockInfo(name=name, content=[line,])
commons.append(c)
return Commons(commons)
def find_names(string):
"""
Find the names in expression ie remove (,),+,-,*,/,**,=
"""
m = re.findall('[A-Z][A-Z0-9_]*', string, re.I)
return set(m)
def find_type(string):
"""
return a Variables object if string is a type declaration
"""
# get out if string is a comment
if string.upper().startswith('C'):
return None
pat = "^\s+"
pat += "(BYTE|"
pat += "CHARACTER(?:\*[0-9]+)?|CHARACTER\*\(\*\)|"
pat += "COMPLEX(?:\*(?:8|16|32))?|"
pat += "DOUBLE COMPLEX|"
pat += "DOUBLE PRECISION|"
pat += "INTEGER(?:\*(?:2|4|8))?|"
pat += "LOGICAL(?:\*(?:1|2|4|8))?|"
pat += "REAL(?:\*(?:4|8|16))?|"
pat += "AUTOMATIC|"
pat += "STATIC)"
pat += "\s+(.*)$"
m = re.match(pat, string, re.I)
if m is not None:
if re.search("IMPLICIT", string):
return None
else:
var_list = find_variables(m.groups()[1])
for var in var_list:
var.type = m.groups()[0]
return var_list
def find_dim(string):
"""
return a Variables object if string is a dimension declaration
"""
pat = "^\s+DIMENSION\s+(.*)$"
m = re.match(pat, string, re.I)
if m is not None:
var_list = find_variables(m.groups()[0])
return var_list
def write_modules(infile):
fi = FileInfo(filename=infile)
# Get all the common blocks defined in the source file
all_commons = []
for sp in fi.subprograms:
for c in sp.commons:
if c.name not in [_.name for _ in all_commons]:
all_commons.append(c)
all_commons = Commons(all_commons)
# a function to create a module fortran code from a CommonBlockInfo object
def common2module(cbi):
variables = cbi.variables
alloc_args = set()
module_name = cbi.name.upper() + "_MOD"
for variable in variables:
if variable.dimension is not None:
dim_vars = find_names(str(variable.dimension))
alloc_args.update(dim_vars)
s = f"MODULE {module_name}\n"
s += " IMPLICIT NONE\n"
# for each variable whose type is defined explicitely
for variable in variables:
s += f" {variable.type}"
dimension = variable.dimension
if dimension is not None:
if dimension.rank > 0:
s += f", ALLOCATABLE, DIMENSION(:" + ",:" * (variable.dimension.rank-1) + ")"
s += f" :: {variable.name}\n"
s += "CONTAINS\n"
#s += f" SUBROUTINE ALLOC_{cbi.name.upper()}({','.join(alloc_args)})\n"
#s += " IMPLICIT INTEGER (A-Z)\n"
s += f" SUBROUTINE ALLOC_{cbi.name.upper()}()\n"
s += f" USE DIM_MOD\n"
# for each variable with a defined dimension
for variable in variables:
if variable.dimension is not None:
s += f" IF (ALLOCATED({variable.name})) THEN\n"
s += f" DEALLOCATE({variable.name})\n"
s += f" ENDIF\n"
s += f" ALLOCATE({variable})\n"
s += f" END SUBROUTINE ALLOC_{cbi.name.upper()}\n"
s += f"END MODULE {module_name}\n"
# indentation
s = textwrap.indent(s, prefix=" ")
# split in too long
content = s.split('\n')
s = content2str(content)
return s
# write the modules.f file
with open("modules.f", "w") as fd:
for c in all_commons:
fd.write("C" + "="*71 + "\n")
s = common2module(c)
fd.write(s)
fd.write("\n"*2)
if __name__ == "__main__":
#infile = 'inv_mat_ms2_la.f'
infile = sys.argv[1]
# write the modules.f file
write_modules(infile)
#exit()
fi = FileInfo(filename=infile)
# Get all the common blocks defined in the source file
all_commons = []
for sp in fi.subprograms:
for c in sp.commons:
if c.name not in [_.name for _ in all_commons]:
all_commons.append(c)
all_commons = Commons(all_commons)
content = fi.content
# write the allocation.f file
dim_vars = [
"NATP_M",
"NATCLU_M_",
"NAT_EQ_M",
"N_CL_L_M",
"NE_M",
"NL_M",
"LI_M",
"NEMET_M",
"NO_ST_M",
"NDIF_M",
"NSO_M",
"NTEMP_M",
"NODES_EX_M",
"NSPIN_M",
"NTH_M",
"NPH_M",
"NDIM_M",
"N_TILT_M",
"N_ORD_M",
"NPATH_M",
"NGR_M"]
s = f"SUBROUTINE ALLOCATION({', '.join([v+'_' for v in dim_vars])})\n"
s += f" USE DIM_MOD\n"
s += f" IMPLICIT INTEGER (A-Z)\n"
for v in dim_vars:
s += f" {v} = {v+'_'}\n"
s += f" CALL INIT_DIM()\n"
s += f"END SUBROUTINE ALLOCATION\n"
# indentation
s = textwrap.indent(s, prefix=" ")
# split in too long
s = content2str(s.split('\n'))
with open("allocation.f", "w") as fd:
fd.write(s)
#exit()
# remove type definitions for variables that are in commons
for sp in fi.subprograms:
# get the list all all variables in all commons in this subprogram
vlist = []
for c in sp.commons:
for v in c.variables:
vlist.append(v.name)
for iline, line in enumerate(sp.content):
print(f"{sp.l0+iline:05d}: {line}")
newline = ''
# comment INCLUDE statement
if re.search('INCLUDE.*spec.inc', line, re.I):
newline = "C" + line
# replace the commons by USE statements
m = re.match("^.*COMMON\s+/(.*)/.*$", line, re.I)
if m is not None:
cmn_name = m.groups()[0]
# Here we test if the common variables are the same than the module
cmn_variables = sp.commons[cmn_name].variables
s = f" USE {cmn_name.upper()}_MOD "
modifications = []
for i, v in enumerate(cmn_variables):
original = all_commons[cmn_name].variables[i].name
if v.name != original:
modifications.append(f"{v.name} => {original}")
s += ", ".join(modifications)
newline = s
# Remove type declaration for variables that are now in modules
allv = find_type(line)
newallv = []
line_ = line
if allv is not None:
# Here the line is a declaration statement
# remove every variables that are also in vlist
for v in allv:
if v.name not in vlist:
# keep this variable in the list
newallv.append(str(v))
# if there is no change
if len(allv) == len(newallv):
line_ = ""
# if no more variables are defined, remove the line
elif len(newallv) == 0:
line_ = "C"
else:
line_ = " " + v.type + " " + ",".join(newallv)
newline = line_
# Remove dimension declaration for variables that are now in modules
allv = find_dim(line)
newallv = []
line_ = line
if allv is not None:
# Here the line is a dimension statement
# remove every variables that are also in vlist
for v in allv:
if v.name not in vlist:
# keep this variable in the list
#if v.dimension is not None:
newallv.append(str(v))
# if there is no change
if len(allv) == len(newallv):
line_ = ""
# if no more variables are defined, remove the line
elif len(newallv) == 0:
line_ = "C"
else:
line_ = " DIMENSION " + ",".join(newallv)
newline = line_
if newline != '':
print(sp.l0, iline, sp.l0+iline)
content[sp.l0+iline] = newline
print(f">>> {newline}")
# rewrite the file
with open(infile + ".new", "w") as fd:
fd.write(content2str(content))

785
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/treat_led.f Normal file
View File

@ -0,0 +1,785 @@
C
C=======================================================================
C
SUBROUTINE TREAT_PHD(ISOM,NFICHLEC,JFICH,NP)
C
C This routine sums up the calculations corresponding to different
C absorbers or different planes when this has to be done
C (parameter ISOM in the input data file).
C
C Last modified : 24 Jan 2013
C
C INCLUDE 'spec.inc'
USE DIM_MOD
USE OUTUNITS_MOD
USE TYPEXP_MOD, DUMMY => SPECTRO
USE VALIN_MOD
USE VALFIN_MOD
C
PARAMETER(N_HEAD=5000,N_FILES=1000)
C
CHARACTER*3 SPECTRO
C
CHARACTER*13 OUTDATA
CHARACTER*72 HEAD(N_HEAD,N_FILES)
C
REAL TAB(NDIM_M,4)
REAL ECIN(NE_M),DTHETA(NTH_M),DPHI(NPH_M)
C
C
DATA JVOL,JTOT/0,-1/
C
REWIND IUO2
C
C Reading and storing the headers:
C
NHEAD=0
DO JLINE=1,N_HEAD
READ(IUO2,888) HEAD(JLINE,JFICH)
NHEAD=NHEAD+1
IF(HEAD(JLINE,JFICH)(1:6).EQ.' ') GOTO 333
ENDDO
C
333 CONTINUE
C
READ(IUO2,15) SPECTRO,OUTDATA
READ(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE,
&IPH_1,I_EXT
C
IF(I_EXT.EQ.2) THEN
IPH_1=0
ENDIF
C
IF(ISOM.EQ.0) THEN
C
C........ ISOM = 0 : case of independent input files .................
C
READ(IUO2,1) NPLAN,NEMET,NTHETA,NPHI,NE
C
IF(IPH_1.EQ.1) THEN
N_FIXED=NPHI
FIX0=PHI0
FIX1=PHI1
N_SCAN=NTHETA
ELSE
N_FIXED=NTHETA
FIX0=THETA0
FIX1=THETA1
IF(STEREO.EQ.'YES') THEN
NPHI=INT((PHI1-PHI0)*FLOAT(NTHETA-1)/(THETA1-THETA0)+
& 0.0001)+1
IF(NTHETA*NPHI.GT.NPH_M) GOTO 37
ENDIF
N_SCAN=NPHI
ENDIF
C
IF(I_EXT.EQ.-1) THEN
N_SCAN=2*N_SCAN
ENDIF
C
IF((I_EXT.EQ.0).OR.(I_EXT.EQ.1)) THEN
NDP=NEMET*NTHETA*NPHI*NE
ELSEIF(I_EXT.EQ.-1) THEN
NDP=NEMET*NTHETA*NPHI*NE*2
ELSEIF(I_EXT.EQ.2) THEN
NDP=NEMET*NTHETA*NE
N_FIXED=NTHETA
N_SCAN=NPHI
IF((N_FIXED.GT.NTH_M).OR.(N_FIXED.GT.NPH_M)) GOTO 35
ENDIF
C
NTT=NPLAN*NDP
IF(NTT.GT.NDIM_M) GOTO 5
C
DO JPLAN=1,NPLAN
DO JEMET=1,NEMET
DO JE=1,NE
C
DO J_FIXED=1,N_FIXED
IF(N_FIXED.GT.1) THEN
XINCRF=FLOAT(J_FIXED-1)*(FIX1-FIX0)/FLOAT(
& N_FIXED-1)
ELSEIF(N_FIXED.EQ.1) THEN
XINCRF=0.
ENDIF
IF(IPH_1.EQ.1) THEN
JPHI=J_FIXED
ELSE
THETA=THETA0+XINCRF
JTHETA=J_FIXED
IF((ABS(THETA).GT.90.).AND.(I_EXT.NE.2)) GOTO 11
ENDIF
IF(STEREO.EQ.' NO') THEN
N_SCAN_R=N_SCAN
ELSE
RTHETA=THETA*0.017453
FIX_STEP=(FIX1-FIX0)/FLOAT(N_FIXED-1)
N_SCAN_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.
& 0001)+1
ENDIF
C
DO J_SCAN=1,N_SCAN_R
IF(IPH_1.EQ.1) THEN
JTHETA=J_SCAN
ELSE
JPHI=J_SCAN
ENDIF
C
JLIN=(JPLAN-1)*NDP + (JEMET-1)*NE*N_FIXED*N_SCAN +
& (JE-1)*N_FIXED*N_SCAN +(JTHETA-1)*NPHI + JPHI
C
IF(I_EXT.LE.0) THEN
IF(STEREO.EQ.' NO') THEN
JPHI2=JPHI
ELSE
JPHI2=(JTHETA-1)*NPHI+JPHI
ENDIF
ELSE
JPHI2=JTHETA
ENDIF
C
READ(IUO2,2) JPL
IF(JPLAN.EQ.JPL) THEN
BACKSPACE IUO2
IF(IDICHR.EQ.0) THEN
READ(IUO2,2) JPL,JEM,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),TAB(JLIN,1),TAB(JLIN,2)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
READ(IUO2,25) TAB(JLIN2,1),TAB(JLIN2,2)
ENDIF
ELSE
READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),TAB(JLIN,1),TAB(JLIN,2),TAB(JLIN,3),
& TAB(JLIN,4)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),
& DPHI(JPHI2),ECIN(JE),TAB(JLIN2,1),TAB(JLIN2,2),TAB(
& JLIN2,3),TAB(JLIN2,4)
ENDIF
ENDIF
ELSE
BACKSPACE IUO2
DO JL=JLIN,JPLAN*NDP
TAB(JL,1)=0.0
TAB(JL,2)=0.0
TAB(JL,3)=0.0
TAB(JL,4)=0.0
ENDDO
GOTO 10
ENDIF
ENDDO
ENDDO
11 CONTINUE
ENDDO
ENDDO
10 CONTINUE
ENDDO
C
REWIND IUO2
C
C Skipping the NHEAD lines of headers before rewriting:
C
DO JLINE=1,NHEAD
READ(IUO2,888) HEAD(JLINE,JFICH)
ENDDO
C
WRITE(IUO2,15) SPECTRO,OUTDATA
WRITE(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
WRITE(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
C
DO JE=1,NE
DO JTHETA=1,NTHETA
IF(STEREO.EQ.' NO') THEN
NPHI_R=NPHI
ELSE
RTHETA=DTHETA(JTHETA)*0.017453
FIX_STEP=(THETA1-THETA0)/FLOAT(NTHETA-1)
NPHI_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.0001)+1
NPHI=INT((PHI1-PHI0)/FIX_STEP+0.0001)+1
ENDIF
DO JPHI=1,NPHI_R
TOTDIF_1=0.
TOTDIR_1=0.
VOLDIF_1=0.
VOLDIR_1=0.
TOTDIF_2=0.
TOTDIR_2=0.
VOLDIF_2=0.
VOLDIR_2=0.
IF(I_EXT.EQ.-1) THEN
TOTDIF2_1=0.
TOTDIR2_1=0.
VOLDIF2_1=0.
VOLDIR2_1=0.
TOTDIF2_2=0.
TOTDIR2_2=0.
VOLDIF2_2=0.
VOLDIR2_2=0.
ENDIF
C
DO JPLAN=1,NPLAN
C
SF_1=0.
SR_1=0.
SF_2=0.
SR_2=0.
IF(I_EXT.EQ.-1) THEN
SF2_1=0.
SR2_1=0.
SF2_2=0.
SR2_2=0.
ENDIF
C
DO JEMET=1,NEMET
JLIN=(JPLAN-1)*NDP + (JEMET-1)*NE*NTHETA*NPHI + (
& JE-1)*NTHETA*NPHI +(JTHETA-1)*NPHI + JPHI
SF_1=SF_1+TAB(JLIN,2)
SR_1=SR_1+TAB(JLIN,1)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
SF2_1=SF2_1+TAB(JLIN2,2)
SR2_1=SR2_1+TAB(JLIN2,1)
ENDIF
IF(IDICHR.GE.1) THEN
SF_2=SF_2+TAB(JLIN,4)
SR_2=SR_2+TAB(JLIN,3)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
SF2_2=SF2_2+TAB(JLIN2,4)
SR2_2=SR2_2+TAB(JLIN2,3)
ENDIF
ENDIF
ENDDO
IF(I_EXT.LE.0) THEN
IF(STEREO.EQ.' NO') THEN
JPHI2=JPHI
ELSE
JPHI2=(JTHETA-1)*NPHI+JPHI
ENDIF
ELSE
JPHI2=JTHETA
ENDIF
IF(IDICHR.EQ.0) THEN
WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),SR_1,SF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR2_1,SF2_1
ENDIF
ELSE
WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),SR_1,SF_1,SR_2,SF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR2_1,SF2_1,SR2_2,SF2_2
ENDIF
ENDIF
IF(JPLAN.GT.NONVOL(JFICH)) THEN
VOLDIF_1=VOLDIF_1+SF_1
VOLDIR_1=VOLDIR_1+SR_1
IF(I_EXT.EQ.-1) THEN
VOLDIF2_1=VOLDIF2_1+SF2_1
VOLDIR2_1=VOLDIR2_1+SR2_1
ENDIF
IF(IDICHR.GE.1) THEN
VOLDIF_2=VOLDIF_2+SF_2
VOLDIR_2=VOLDIR_1+SR_2
IF(I_EXT.EQ.-1) THEN
VOLDIF2_2=VOLDIF2_2+SF2_2
VOLDIR2_2=VOLDIR2_1+SR2_2
ENDIF
ENDIF
ENDIF
TOTDIF_1=TOTDIF_1+SF_1
TOTDIR_1=TOTDIR_1+SR_1
IF(I_EXT.EQ.-1) THEN
TOTDIF2_1=TOTDIF2_1+SF2_1
TOTDIR2_1=TOTDIR2_1+SR2_1
ENDIF
IF(IDICHR.GE.1) THEN
TOTDIF_2=TOTDIF_2+SF_2
TOTDIR_2=TOTDIR_2+SR_2
IF(I_EXT.EQ.-1) THEN
TOTDIF2_2=TOTDIF2_2+SF2_2
TOTDIR2_2=TOTDIR2_2+SR2_2
ENDIF
ENDIF
ENDDO
IF(IDICHR.EQ.0) THEN
WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(
& JE),VOLDIR_1,VOLDIF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),VOLDIR2_1,VOLDIF2_1
ENDIF
WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(
& JE),TOTDIR_1,TOTDIF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),TOTDIR2_1,TOTDIF2_1
ENDIF
ELSE
WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),VOLDIR_1,VOLDIF_1,VOLDIR_2,VOLDIF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),VOLDIR2_1,VOLDIF2_1,VOLDIR2_2,VOLDIF2_2
ENDIF
WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),TOTDIR_1,TOTDIF_1,TOTDIR_2,TOTDIF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),TOTDIR2_1,TOTDIF2_1,TOTDIR2_2,TOTDIF2_2
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
C
ELSE
C
C........ ISOM not= 0 : multiple input files to be summed up ..........
C
READ(IUO2,7) NTHETA,NPHI,NE
C
IF(IPH_1.EQ.1) THEN
N_FIXED=NPHI
FIX0=PHI0
FIX1=PHI1
N_SCAN=NTHETA
ELSE
N_FIXED=NTHETA
FIX0=THETA0
FIX1=THETA1
IF(STEREO.EQ.'YES') THEN
NPHI=INT((PHI1-PHI0)*FLOAT(NTHETA-1)/(THETA1-THETA0)+
& 0.0001)+1
IF(NTHETA*NPHI.GT.NPH_M) GOTO 37
ENDIF
N_SCAN=NPHI
ENDIF
C
IF(I_EXT.EQ.-1) THEN
N_SCAN=2*N_SCAN
ENDIF
C
IF((I_EXT.EQ.0).OR.(I_EXT.EQ.1)) THEN
NDP=NTHETA*NPHI*NE
ELSEIF(I_EXT.EQ.-1) THEN
NDP=NTHETA*NPHI*NE*2
ELSEIF(I_EXT.EQ.2) THEN
NDP=NTHETA*NE
N_FIXED=NTHETA
N_SCAN=NPHI
IF((N_FIXED.GT.NTH_M).OR.(N_FIXED.GT.NPH_M)) GOTO 35
ENDIF
C
NTT=NFICHLEC*NDP
IF(NTT.GT.NDIM_M) GOTO 5
C
IF(ISOM.EQ.1) THEN
NPLAN=NP
NF=NP
ELSEIF(ISOM.EQ.2) THEN
NEMET=NFICHLEC
NF=NFICHLEC
NPLAN=1
ENDIF
C
DO JF=1,NF
C
C Reading the headers for each file:
C
IF(JF.GT.1) THEN
DO JLINE=1,NHEAD
READ(IUO2,888) HEAD(JLINE,JF)
ENDDO
ENDIF
C
DO JE=1,NE
C
DO J_FIXED=1,N_FIXED
IF(N_FIXED.GT.1) THEN
XINCRF=FLOAT(J_FIXED-1)*(FIX1-FIX0)/FLOAT(
& N_FIXED-1)
ELSEIF(N_FIXED.EQ.1) THEN
XINCRF=0.
ENDIF
IF(IPH_1.EQ.1) THEN
JPHI=J_FIXED
ELSE
THETA=THETA0+XINCRF
JTHETA=J_FIXED
IF((ABS(THETA).GT.90.).AND.(I_EXT.NE.2)) GOTO 12
ENDIF
IF(STEREO.EQ.' NO') THEN
N_SCAN_R=N_SCAN
ELSE
RTHETA=THETA*0.017453
FIX_STEP=(FIX1-FIX0)/FLOAT(N_FIXED-1)
N_SCAN_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.
& 0001)+1
ENDIF
C
DO J_SCAN=1,N_SCAN_R
IF(IPH_1.EQ.1) THEN
JTHETA=J_SCAN
ELSE
JPHI=J_SCAN
ENDIF
C
JLIN=(JF-1)*NDP + (JE-1)*N_FIXED*N_SCAN +(JTHETA-1)
& *NPHI + JPHI
IF(I_EXT.LE.0) THEN
IF(STEREO.EQ.' NO') THEN
JPHI2=JPHI
ELSE
JPHI2=(JTHETA-1)*NPHI+JPHI
ENDIF
ELSE
JPHI2=JTHETA
ENDIF
C
IF(ISOM.EQ.1) THEN
READ(IUO2,2) JPL
IF(JF.EQ.JPL) THEN
BACKSPACE IUO2
IF(IDICHR.EQ.0) THEN
READ(IUO2,2) JPL,JEM,DTHETA(JTHETA),
& DPHI(JPHI2),ECIN(JE),TAB(JLIN,1),TAB(JLIN,2)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
READ(IUO2,25) TAB(JLIN2,1),TAB(
& JLIN2,2)
ENDIF
ELSE
READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),
& DPHI(JPHI2),ECIN(JE),TAB(JLIN,1),TAB(JLIN,2),TAB(
& JLIN,3),TAB(JLIN,4)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
READ(IUO2,22) JPL,JEM,DTHETA(
& JTHETA),DPHI(JPHI2),ECIN(JE),TAB(JLIN2,1),TAB(
& JLIN2,2),TAB(JLIN2,3),TAB(JLIN2,4)
ENDIF
ENDIF
ELSE
BACKSPACE IUO2
DO JLINE=1,NHEAD
BACKSPACE IUO2
ENDDO
DO JL=JLIN,JF*NDP
TAB(JL,1)=0.0
TAB(JL,2)=0.0
TAB(JL,3)=0.0
TAB(JL,4)=0.0
ENDDO
GOTO 13
ENDIF
ELSEIF(ISOM.EQ.2) THEN
IF(IDICHR.EQ.0) THEN
READ(IUO2,2) JPL,JEM,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),TAB(JLIN,1),TAB(JLIN,2)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
READ(IUO2,25) TAB(JLIN2,1),TAB(JLIN2,2)
ENDIF
ELSE
READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),TAB(JLIN,1),TAB(JLIN,2),TAB(JLIN,3),
& TAB(JLIN,4)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
READ(IUO2,22) JPL,JEM,DTHETA(JTHETA),
& DPHI(JPHI2),ECIN(JE),TAB(JLIN2,1),TAB(JLIN2,2),TAB(
& JLIN2,3),TAB(JLIN2,4)
ENDIF
ENDIF
ENDIF
ENDDO
12 CONTINUE
ENDDO
ENDDO
13 CONTINUE
ENDDO
C
REWIND IUO2
C
C Writing the headers:
C
DO JLINE=1,2
WRITE(IUO2,888) HEAD(JLINE,1)
ENDDO
DO JF=1,NFICHLEC
DO JLINE=3,6
WRITE(IUO2,888) HEAD(JLINE,JF)
ENDDO
WRITE(IUO2,888) HEAD(2,JF)
ENDDO
DO JLINE=7,NHEAD
WRITE(IUO2,888) HEAD(JLINE,1)
ENDDO
C
WRITE(IUO2,15) SPECTRO,OUTDATA
WRITE(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
WRITE(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
C
IF(ISOM.EQ.1) THEN
C
DO JE=1,NE
C
DO JTHETA=1,NTHETA
IF(STEREO.EQ.' NO') THEN
NPHI_R=NPHI
ELSE
RTHETA=DTHETA(JTHETA)*0.017453
FIX_STEP=(THETA1-THETA0)/FLOAT(NTHETA-1)
NPHI_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.
& 0001)+1
NPHI=INT((PHI1-PHI0)/FIX_STEP+0.0001)+1
ENDIF
DO JPHI=1,NPHI_R
C
TOTDIF_1=0.
TOTDIR_1=0.
VOLDIF_1=0.
VOLDIR_1=0.
TOTDIF_2=0.
TOTDIR_2=0.
VOLDIF_2=0.
VOLDIR_2=0.
IF(I_EXT.EQ.-1) THEN
TOTDIF2_1=0.
TOTDIR2_1=0.
VOLDIF2_1=0.
VOLDIR2_1=0.
TOTDIF2_2=0.
TOTDIR2_2=0.
VOLDIF2_2=0.
VOLDIR2_2=0.
ENDIF
C
DO JPLAN=1,NPLAN
JF=JPLAN
C
JLIN=(JF-1)*NDP + (JE-1)*NTHETA*NPHI +(JTHETA-1)*
& NPHI + JPHI
C
SR_1=TAB(JLIN,1)
SF_1=TAB(JLIN,2)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
SF2_1=TAB(JLIN2,2)
SR2_1=TAB(JLIN2,1)
ENDIF
IF(I_EXT.LE.0) THEN
IF(STEREO.EQ.' NO') THEN
JPHI2=JPHI
ELSE
JPHI2=(JTHETA-1)*NPHI+JPHI
ENDIF
ELSE
JPHI2=JTHETA
ENDIF
IF(IDICHR.EQ.0) THEN
WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR_1,SF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR2_1,SF2_1
ENDIF
ELSE
SR_2=TAB(JLIN,3)
SF_2=TAB(JLIN,4)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
SF2_2=TAB(JLIN2,4)
SR2_2=TAB(JLIN2,3)
ENDIF
WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR_1,SF_1,SR_2,SF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR2_1,SF2_1,SR2_2,SF2_2
ENDIF
ENDIF
IF(NONVOL(JPLAN).EQ.0) THEN
VOLDIF_1=VOLDIF_1+SF_1
VOLDIR_1=VOLDIR_1+SR_1
IF(I_EXT.EQ.-1) THEN
VOLDIF2_1=VOLDIF2_1+SF2_1
VOLDIR2_1=VOLDIR2_1+SR2_1
ENDIF
IF(IDICHR.GE.1) THEN
VOLDIF_2=VOLDIF_2+SF_2
VOLDIR_2=VOLDIR_2+SR_2
IF(I_EXT.EQ.-1) THEN
VOLDIF2_2=VOLDIF2_2+SF2_2
VOLDIR2_2=VOLDIR2_1+SR2_2
ENDIF
ENDIF
ENDIF
TOTDIF_1=TOTDIF_1+SF_1
TOTDIR_1=TOTDIR_1+SR_1
IF(I_EXT.EQ.-1) THEN
TOTDIF2_1=TOTDIF2_1+SF2_1
TOTDIR2_1=TOTDIR2_1+SR2_1
ENDIF
IF(IDICHR.GE.1) THEN
TOTDIF_2=TOTDIF_2+SF_2
TOTDIR_2=TOTDIR_2+SR_2
IF(I_EXT.EQ.-1) THEN
TOTDIF2_2=TOTDIF2_2+SF2_2
TOTDIR2_2=TOTDIR2_2+SR2_2
ENDIF
ENDIF
ENDDO
C
IF(IDICHR.EQ.0) THEN
WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),VOLDIR_1,VOLDIF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),VOLDIR2_1,VOLDIF2_1
ENDIF
WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),TOTDIR_1,TOTDIF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),TOTDIR2_1,TOTDIF2_1
ENDIF
ELSE
WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),VOLDIR_1,VOLDIF_1,VOLDIR_2,VOLDIF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),VOLDIR2_1,VOLDIF2_1,VOLDIR2_2,
& VOLDIF2_2
ENDIF
WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),TOTDIR_1,TOTDIF_1,TOTDIR_2,TOTDIF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),TOTDIR2_1,TOTDIF2_1,TOTDIR2_2,
& TOTDIF2_2
ENDIF
ENDIF
C
ENDDO
ENDDO
ENDDO
ELSEIF(ISOM.EQ.2) THEN
DO JE=1,NE
C
DO JTHETA=1,NTHETA
IF(STEREO.EQ.' NO') THEN
NPHI_R=NPHI
ELSE
RTHETA=DTHETA(JTHETA)*0.017453
FIX_STEP=(THETA1-THETA0)/FLOAT(NTHETA-1)
NPHI_R=INT((PHI1-PHI0)*SIN(RTHETA)/FIX_STEP+0.
& 0001)+1
NPHI=INT((PHI1-PHI0)/FIX_STEP+0.0001)+1
ENDIF
DO JPHI=1,NPHI_R
C
SF_1=0.
SR_1=0.
SF_2=0.
SR_2=0.
IF(I_EXT.EQ.-1) THEN
SF2_1=0.
SR2_1=0.
SF2_2=0.
SR2_2=0.
ENDIF
C
DO JEMET=1,NEMET
JF=JEMET
C
JLIN=(JF-1)*NDP + (JE-1)*NTHETA*NPHI +(JTHETA-
& 1)*NPHI + JPHI
C
SF_1=SF_1+TAB(JLIN,2)
SR_1=SR_1+TAB(JLIN,1)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
SF2_1=SF2_1+TAB(JLIN2,2)
SR2_1=SR2_1+TAB(JLIN2,1)
ENDIF
IF(IDICHR.GE.1) THEN
SF_2=SF_2+TAB(JLIN,4)
SR_2=SR_2+TAB(JLIN,3)
IF(I_EXT.EQ.-1) THEN
JLIN2=NTT+JLIN
SF2_2=SF2_2+TAB(JLIN2,4)
SR2_2=SR2_2+TAB(JLIN2,3)
ENDIF
ENDIF
ENDDO
IF(I_EXT.LE.0) THEN
IF(STEREO.EQ.' NO') THEN
JPHI2=JPHI
ELSE
JPHI2=(JTHETA-1)*NPHI+JPHI
ENDIF
ELSE
JPHI2=JTHETA
ENDIF
IF(IDICHR.EQ.0) THEN
WRITE(IUO2,3) JPL,DTHETA(JTHETA),DPHI(JPHI2),
& ECIN(JE),SR_1,SF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR2_1,SF2_1
ENDIF
ELSE
WRITE(IUO2,23) JPL,DTHETA(JTHETA),DPHI(JPHI2)
& ,ECIN(JE),SR_1,SF_1,SR_2,SF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(
& JPHI2),ECIN(JE),SR2_1,SF2_1,SR2_2,SF2_2
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
ENDIF
C
GOTO 6
C
5 WRITE(IUO1,4)
STOP
35 WRITE(IUO1,36) N_FIXED
STOP
37 WRITE(IUO1,38) NTHETA*NPHI
STOP
C
1 FORMAT(2X,I3,2X,I2,2X,I4,2X,I4,2X,I4)
2 FORMAT(2X,I3,2X,I2,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
3 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
4 FORMAT(//,8X,'<<<<<<<<<< DIMENSION OF THE ARRAYS TOO SMALL ',
&'IN THE TREAT_PHD SUBROUTINE - INCREASE NDIM_M ','>>>>>>>>>>')
7 FORMAT(I4,2X,I4,2X,I4)
8 FORMAT(I4,2X,I4,2X,I4,2X,I3,2X,I1)
9 FORMAT(9(2X,I1),2X,I2)
15 FORMAT(2X,A3,11X,A13)
22 FORMAT(2X,I3,2X,I2,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,
&E12.6,2X,E12.6)
23 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,
&2X,E12.6)
25 FORMAT(37X,E12.6,2X,E12.6)
36 FORMAT(//,4X,'<<<<<<<<<< DIMENSION OF NTH_M OR NPH_M TOO SMALL
&','IN THE INCLUDE FILE >>>>>>>>>>',/,4X,'<<<<<<<<<<
& SHOULD BE AT LEAST ',I6,' >>>>>>>>>>')
38 FORMAT(//,8X,'<<<<<<<<<< DIMENSION OF NPH_M TOO SMALL ','IN THE
&INCLUDE FILE >>>>>>>>>>',/,8X,'<<<<<<<<<< SHOULD BE
&AT LEAST ',I6,' >>>>>>>>>>')
888 FORMAT(A72)
C
6 RETURN
C
END

335
src/msspec/spec/fortran/led_mi_noso_nosp_nosym/weight_sum.f Normal file
View File

@ -0,0 +1,335 @@
C
C=======================================================================
C
SUBROUTINE WEIGHT_SUM(ISOM,I_EXT,I_EXT_A,JEL)
C
C This subroutine performs a weighted sum of the results
C corresponding to different directions of the detector.
C The directions and weights are read from an external input file
C
C JEL is the electron undetected (i.e. for which the outgoing
C directions are integrated over the unit sphere). It is always
C 1 for one electron spectroscopies (PHD). For APECS, It can be
C 1 (photoelectron) or 2 (Auger electron) or even 0 (no electron
C detected)
C
C Last modified : 31 Jan 2007
C
USE DIM_MOD
USE INFILES_MOD
USE INUNITS_MOD
USE OUTUNITS_MOD
C
C
PARAMETER(N_MAX=5810,NPM=20)
C
REAL*4 W(N_MAX),W_A(N_MAX),ECIN(NE_M)
REAL*4 DTHETA(N_MAX),DPHI(N_MAX),DTHETAA(N_MAX),DPHIA(N_MAX)
REAL*4 SR_1,SF_1,SR_2,SF_2
REAL*4 SUMR_1(NPM,NE_M,N_MAX),SUMR_2(NPM,NE_M,N_MAX)
REAL*4 SUMF_1(NPM,NE_M,N_MAX),SUMF_2(NPM,NE_M,N_MAX)
C
CHARACTER*3 SPECTRO,SPECTRO2
CHARACTER*5 LIKE
CHARACTER*13 OUTDATA
C
C
C
C
DATA JVOL,JTOT/0,-1/
DATA LIKE /'-like'/
C
REWIND IUO2
C
READ(IUO2,15) SPECTRO,OUTDATA
IF(SPECTRO.NE.'APC') THEN
READ(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
READ(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
SPECTRO2='XAS'
ELSE
READ(IUO2,9) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
READ(IUO2,9) ISPIN_A,IDICHR_A,I_SO_A,ISFLIP_A,ICHKDIR_A,IPHI_A,I
&THETA_A,IE_A
READ(IUO2,8) NPHI,NTHETA,NE,NPLAN,ISOM
READ(IUO2,8) NPHI_A,NTHETA_A
IF(JEL.EQ.1) THEN
SPECTRO2='AED'
ELSEIF(JEL.EQ.2) THEN
SPECTRO2='PHD'
ELSEIF(JEL.EQ.0) THEN
SPECTRO2='XAS'
ENDIF
ENDIF
C
IF(NPLAN.GT.NPM) THEN
WRITE(IUO1,4) NPLAN+2
STOP
ENDIF
C
C Reading the number of angular points
C
IF(SPECTRO.NE.'APC') THEN
OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
READ(IUI6,1) N_POINTS
READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
N_POINTS_A=1
ELSE
IF(JEL.EQ.1) THEN
OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
READ(IUI6,1) N_POINTS
READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
IF(I_EXT_A.EQ.0) THEN
N_POINTS_A=NTHETA_A*NPHI_A
ELSE
OPEN(UNIT=IUI9, FILE=INFILE9, STATUS='OLD')
READ(IUI9,1) N_POINTS_A
READ(IUI9,5) I_DIM,N_DUM1,N_DUM2
ENDIF
NTHETA0=NTHETA_A
NPHI0=NPHI_A
ELSEIF(JEL.EQ.2) THEN
OPEN(UNIT=IUI9, FILE=INFILE9, STATUS='OLD')
READ(IUI9,1) N_POINTS_A
READ(IUI9,5) I_DIM,N_DUM1,N_DUM2
IF(I_EXT.EQ.0) THEN
N_POINTS=NTHETA*NPHI
ELSE
OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
READ(IUI6,1) N_POINTS
READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
ENDIF
NTHETA0=NTHETA
NPHI0=NPHI
ELSEIF(JEL.EQ.0) THEN
OPEN(UNIT=IUI6, FILE=INFILE6, STATUS='OLD')
OPEN(UNIT=IUI9, FILE=INFILE9, STATUS='OLD')
READ(IUI6,1) N_POINTS
READ(IUI9,1) N_POINTS_A
READ(IUI6,5) I_DIM,N_DUM1,N_DUM2
READ(IUI9,5) I_DIM,N_DUM1,N_DUM2
ENDIF
ENDIF
C
IF(SPECTRO.NE.'APC') THEN
NANGLE=1
ELSE
IF(JEL.EQ.1) THEN
NANGLE=N_POINTS_A
ELSEIF(JEL.EQ.2) THEN
NANGLE=N_POINTS
ELSEIF(JEL.EQ.0) THEN
NANGLE=1
ENDIF
ENDIF
C
C Initialization of the arrays
C
DO JE=1,NE
DO JANGLE=1,NANGLE
DO JPLAN=1,NPLAN+2
SUMR_1(JPLAN,JE,JANGLE)=0.
SUMF_1(JPLAN,JE,JANGLE)=0.
IF(IDICHR.GT.0) THEN
SUMR_2(JPLAN,JE,JANGLE)=0.
SUMF_2(JPLAN,JE,JANGLE)=0.
ENDIF
ENDDO
ENDDO
ENDDO
C
C Reading of the data to be angle integrated
C
DO JE=1,NE
C
DO JANGLE=1,N_POINTS
IF(I_EXT.NE.0) READ(IUI6,2) TH,PH,W(JANGLE)
DO JANGLE_A=1,N_POINTS_A
IF((I_EXT_A.NE.0).AND.(JANGLE.EQ.1)) THEN
READ(IUI9,2) THA,PHA,W_A(JANGLE_A)
ENDIF
C
DO JPLAN=1,NPLAN+2
C
IF(IDICHR.EQ.0) THEN
IF(SPECTRO.NE.'APC') THEN
READ(IUO2,3) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE)
&,SR_1,SF_1
ELSE
READ(IUO2,13) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE
&),DTHETAA(JANGLE_A),DPHIA(JANGLE_A),SR_1,SF_1
ENDIF
ELSE
IF(SPECTRO.NE.'APC') THEN
READ(IUO2,23) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE
&),SR_1,SF_1,SR_2,SF_2
ELSE
READ(IUO2,24) JDUM,DTHETA(JANGLE),DPHI(JANGLE),ECIN(JE
&),DTHETAA(JANGLE_A),DPHIA(JANGLE_A),SR_1,SF_1,SR_2,SF_2
ENDIF
ENDIF
C
IF(JEL.EQ.1) THEN
SUMR_1(JPLAN,JE,JANGLE_A)=SUMR_1(JPLAN,JE,JANGLE_A)+SR_1
&*W(JANGLE)
SUMF_1(JPLAN,JE,JANGLE_A)=SUMF_1(JPLAN,JE,JANGLE_A)+SF_1
&*W(JANGLE)
ELSEIF(JEL.EQ.2) THEN
SUMR_1(JPLAN,JE,JANGLE)=SUMR_1(JPLAN,JE,JANGLE)+SR_1*W_A
&(JANGLE_A)
SUMF_1(JPLAN,JE,JANGLE)=SUMF_1(JPLAN,JE,JANGLE)+SF_1*W_A
&(JANGLE_A)
ELSEIF(JEL.EQ.0) THEN
SUMR_1(JPLAN,JE,1)=SUMR_1(JPLAN,JE,1)+SR_1*W(JANGLE)*W_A
&(JANGLE_A)
SUMF_1(JPLAN,JE,1)=SUMF_1(JPLAN,JE,1)+SF_1*W(JANGLE)*W_A
&(JANGLE_A)
ENDIF
IF(IDICHR.GT.0) THEN
IF(JEL.EQ.1) THEN
SUMR_2(JPLAN,JE,JANGLE_A)=SUMR_2(JPLAN,JE,JANGLE_A)+SR
&_2*W(JANGLE)
SUMF_2(JPLAN,JE,JANGLE_A)=SUMF_2(JPLAN,JE,JANGLE_A)+SF
&_2*W(JANGLE)
ELSEIF(JEL.EQ.2) THEN
SUMR_2(JPLAN,JE,JANGLE)=SUMR_2(JPLAN,JE,JANGLE)+SR_2*W
&_A(JANGLE_A)
SUMF_2(JPLAN,JE,JANGLE)=SUMF_2(JPLAN,JE,JANGLE)+SF_2*W
&_A(JANGLE_A)
ELSEIF(JEL.EQ.0) THEN
SUMR_2(JPLAN,JE,1)=SUMR_2(JPLAN,JE,1)+SR_2*W(JANGLE)*W
&_A(JANGLE_A)
SUMF_2(JPLAN,JE,1)=SUMF_2(JPLAN,JE,1)+SF_2*W(JANGLE)*W
&_A(JANGLE_A)
ENDIF
ENDIF
C
ENDDO
C
ENDDO
IF(I_EXT_A.NE.0) THEN
REWIND IUI9
READ(IUI9,1) NDUM
READ(IUI9,1) NDUM
ENDIF
ENDDO
C
IF(I_EXT.NE.0) THEN
REWIND IUI6
READ(IUI6,1) NDUM
READ(IUI6,1) NDUM
ENDIF
ENDDO
C
CLOSE(IUI6)
CLOSE(IUI9)
REWIND IUO2
C
WRITE(IUO2,16) SPECTRO2,LIKE,SPECTRO,OUTDATA
IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
WRITE(IUO2,19) ISPIN,IDICHR,I_SO,ISFLIP
WRITE(IUO2,18) NE,NPLAN,ISOM
ELSEIF(JEL.EQ.1) THEN
WRITE(IUO2,20) ISPIN_A,IDICHR_A,I_SO_A,ISFLIP_A,ICHKDIR_A,IPHI_A
&,ITHETA_A,IE_A
WRITE(IUO2,21) NPHI0,NTHETA0,NE,NPLAN,ISOM
ELSEIF(JEL.EQ.2) THEN
WRITE(IUO2,20) ISPIN,IDICHR,I_SO,ISFLIP,ICHKDIR,IPHI,ITHETA,IE
WRITE(IUO2,21) NPHI0,NTHETA0,NE,NPLAN,ISOM
ENDIF
C
DO JE=1,NE
DO JANGLE=1,NANGLE
IF(SPECTRO.EQ.'APC') THEN
IF(JEL.EQ.1) THEN
THETA=DTHETAA(JANGLE)
PHI=DPHIA(JANGLE)
ELSEIF(JEL.EQ.2) THEN
THETA=DTHETA(JANGLE)
PHI=DPHI(JANGLE)
ENDIF
ENDIF
C
DO JPLAN=1,NPLAN
IF(IDICHR.EQ.0) THEN
IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
WRITE(IUO2,33) JPLAN,ECIN(JE),SUMR_1(JPLAN,JE,JANGLE),SU
&MF_1(JPLAN,JE,JANGLE)
ELSE
WRITE(IUO2,34) JPLAN,THETA,PHI,ECIN(JE),SUMR_1(JPLAN,JE,
&JANGLE),SUMF_1(JPLAN,JE,JANGLE)
ENDIF
ELSE
IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
WRITE(IUO2,43) JPLAN,ECIN(JE),SUMR_1(JPLAN,JE,JANGLE),SU
&MF_1(JPLAN,JE,JANGLE),SUMR_2(JPLAN,JE,JANGLE),SUMF_2(JPLAN,JE,JANG
&LE)
ELSE
WRITE(IUO2,44) JPLAN,THETA,PHI,ECIN(JE),SUMR_1(JPLAN,JE,
&JANGLE),SUMF_1(JPLAN,JE,JANGLE),SUMR_2(JPLAN,JE,JANGLE),SUMF_2(JPL
&AN,JE,JANGLE)
ENDIF
ENDIF
ENDDO
C
IF(IDICHR.EQ.0) THEN
IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
WRITE(IUO2,33) JVOL,ECIN(JE),SUMR_1(NPLAN+1,JE,JANGLE),SUM
&F_1(NPLAN+1,JE,JANGLE)
WRITE(IUO2,33) JTOT,ECIN(JE),SUMR_1(NPLAN+2,JE,JANGLE),SUM
&F_1(NPLAN+2,JE,JANGLE)
ELSE
WRITE(IUO2,34) JVOL,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+1,JE,J
&ANGLE),SUMF_1(NPLAN+1,JE,JANGLE)
WRITE(IUO2,34) JTOT,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+2,JE,J
&ANGLE),SUMF_1(NPLAN+2,JE,JANGLE)
ENDIF
ELSE
IF((SPECTRO.NE.'APC').OR.(JEL.EQ.0)) THEN
WRITE(IUO2,43) JVOL,ECIN(JE),SUMR_1(NPLAN+1,JE,JANGLE),SUM
&F_1(NPLAN+1,JE,JANGLE),SUMR_2(NPLAN+1,JE,JANGLE),SUMF_2(NPLAN+1,JE
&,JANGLE)
WRITE(IUO2,43) JTOT,ECIN(JE),SUMR_1(NPLAN+2,JE,JANGLE),SUM
&F_1(NPLAN+2,JE,JANGLE),SUMR_2(NPLAN+2,JE,JANGLE),SUMF_2(NPLAN+2,JE
&,JANGLE)
ELSE
WRITE(IUO2,44) JVOL,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+1,JE,J
&ANGLE),SUMF_1(NPLAN+1,JE,JANGLE),SUMR_2(NPLAN+1,JE,JANGLE),SUMF_2(
&NPLAN+1,JE,JANGLE)
WRITE(IUO2,44) JTOT,THETA,PHI,ECIN(JE),SUMR_1(NPLAN+2,JE,J
&ANGLE),SUMF_1(NPLAN+2,JE,JANGLE),SUMR_2(NPLAN+2,JE,JANGLE),SUMF_2(
&NPLAN+2,JE,JANGLE)
ENDIF
ENDIF
C
ENDDO
ENDDO
C
1 FORMAT(13X,I4)
2 FORMAT(15X,F8.3,3X,F8.3,3X,E12.6)
3 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
4 FORMAT(//,8X,'<<<<<<<<<< DIMENSION OF THE ARRAYS TOO SMALL ','IN
&THE WEIGHT_SUM SUBROUTINE - INCREASE NPM TO ',I3,'>>>>>>>>>>')
5 FORMAT(6X,I1,1X,I3,3X,I3)
8 FORMAT(I4,2X,I4,2X,I4,2X,I3,2X,I1)
9 FORMAT(9(2X,I1),2X,I2)
13 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,F6.2,2X,F6.2,2X,E12.6,2X,E
&12.6)
15 FORMAT(2X,A3,11X,A13)
16 FORMAT(2X,A3,A5,1X,A3,2X,A13)
18 FORMAT(I4,2X,I3,2X,I1)
19 FORMAT(4(2X,I1))
20 FORMAT(8(2X,I1))
21 FORMAT(I4,2X,I4,2X,I4,2X,I3,2X,I1)
23 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X
&,E12.6)
24 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,F6.2,2X,F6.2,2X,E12.6,2X,E
&12.6,2X,E12.6,2X,E12.6)
33 FORMAT(2X,I3,2X,F8.2,2X,E12.6,2X,E12.6)
34 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6)
43 FORMAT(2X,I3,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X,E12.6)
44 FORMAT(2X,I3,2X,F6.2,2X,F6.2,2X,F8.2,2X,E12.6,2X,E12.6,2X,E12.6,2X
&,E12.6)
C
RETURN
C
END