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

Add shared object for PED by Matrix inversion. 

The PED spectroscopy by the matrix inversion algorithm is now
included.

This commit also includes:
 * some changes in default values in parameters.py for eigval_ipwm, eigval_method
   in SpecParameters and phi in SourceParameters
 * a modification in the cluster viewer: the mouse wheel allows to change the
   transparency and see the emitter
 * in utils.py, the hemispherical_cluster function has a new keyword that allows
   to create cylindrical clusters
This commit is contained in:
Sylvain Tricot 2020-02-21 15:18:53 +01:00
parent 27c772004c
commit 5a817ab97d
24 changed files with 9803 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/msspec/calculator.py
View File

@ -62,7 +62,7 @@ from msspec.calcio import PhagenIO, SpecIO
from msspec.phagen.fortran.libphagen import main as do_phagen from msspec.phagen.fortran.libphagen import main as do_phagen
from msspec.spec.fortran import phd_se_noso_nosp_nosym from msspec.spec.fortran import phd_se_noso_nosp_nosym, phd_mi_noso_nosp_nosym
from msspec.spec.fortran import eig_mi from msspec.spec.fortran import eig_mi
from msspec.spec.fortran import eig_pw from msspec.spec.fortran import eig_pw
@ -300,7 +300,7 @@ class _MSCALCULATOR(Calculator):
os.chdir(self.init_folder) os.chdir(self.init_folder)
def run_spec(self): def run_spec(self, malloc={}):
def get_li(level): def get_li(level):
orbitals = 'spdfghi' orbitals = 'spdfghi'
m = re.match(r'\d(?P<l>[%s])(\d/2)?' % orbitals, level) m = re.match(r'\d(?P<l>[%s])(\d/2)?' % orbitals, level)
@ -351,6 +351,11 @@ class _MSCALCULATOR(Calculator):
'NPATH_M' : 500, 'NPATH_M' : 500,
'NGR_M' : 10,}) 'NGR_M' : 10,})
# update with provided values
for key, value in malloc.items():
requirements[key] = value
# set some automatic values for memory allocation
for key, value in requirements.items(): for key, value in requirements.items():
setattr(self.spec_malloc_parameters, key, value) setattr(self.spec_malloc_parameters, key, value)
@ -358,6 +363,8 @@ class _MSCALCULATOR(Calculator):
if self.global_parameters.spectroscopy == 'PED': if self.global_parameters.spectroscopy == 'PED':
if self.global_parameters.algorithm == 'expansion': if self.global_parameters.algorithm == 'expansion':
do_spec = phd_se_noso_nosp_nosym.run do_spec = phd_se_noso_nosp_nosym.run
elif self.global_parameters.algorithm == 'inversion':
do_spec = phd_mi_noso_nosp_nosym.run
else: else:
LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not " LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
"an allowed combination.".format(self.global_parameters.spectroscopy, "an allowed combination.".format(self.global_parameters.spectroscopy,
@ -447,7 +454,7 @@ class _MSCALCULATOR(Calculator):
try: try:
# for each set of tl: # for each set of tl:
# 1. get the symbol of the prototipical atom # 1. get the symbol of the prototipical atom
j = np.where(proto_indices == ia+1) j = np.where(proto_indices == ia+1)[0]
symbol = cluster[j][0].symbol symbol = cluster[j][0].symbol
# 2. get the number of max tl allowed # 2. get the number of max tl allowed
ntl = max_tl[symbol] ntl = max_tl[symbol]
@ -561,7 +568,8 @@ class _PED(_MSCALCULATOR):
def _get_scan(self, scan_type='theta', phi=0, def _get_scan(self, scan_type='theta', phi=0,
theta=np.linspace(-70, 70, 141), level=None, theta=np.linspace(-70, 70, 141), level=None,
kinetic_energy=None, data=None): kinetic_energy=None, data=None,
malloc={}):
LOGGER.info("Computting the %s scan...", scan_type) LOGGER.info("Computting the %s scan...", scan_type)
if data: if data:
self.iodata = data self.iodata = data
@ -597,7 +605,7 @@ class _PED(_MSCALCULATOR):
self.spectroscopy_parameters.set_parameter('level', level) self.spectroscopy_parameters.set_parameter('level', level)
self.get_tmatrix() self.get_tmatrix()
self.run_spec() self.run_spec(malloc)
# Now load the data # Now load the data
ndset = len(self.iodata) ndset = len(self.iodata)
@ -858,9 +866,9 @@ class _PED(_MSCALCULATOR):
argument or a new :py:class:`iodata.Data` object. argument or a new :py:class:`iodata.Data` object.
""" """
self.spec_malloc_parameters.NPH_M = 8000
data = self._get_scan(scan_type='theta_phi', level=level, theta=theta, data = self._get_scan(scan_type='theta_phi', level=level, theta=theta,
phi=phi, kinetic_energy=kinetic_energy, data=data) phi=phi, kinetic_energy=kinetic_energy, data=data,
malloc={'NPH_M': 8000})
return data return data

4
src/msspec/iodata.py
View File

@ -573,6 +573,8 @@ class Data(object):
try: try:
del meta_grp['info'] del meta_grp['info']
except: except:
pass
finally:
meta_grp.create_dataset('info', data=np.array((xml_str,)).view('S1')) meta_grp.create_dataset('info', data=np.array((xml_str,)).view('S1'))
self._dirty = False self._dirty = False
LOGGER.info('Data saved in {}'.format(os.path.abspath(filename))) LOGGER.info('Data saved in {}'.format(os.path.abspath(filename)))
@ -998,7 +1000,7 @@ class _DataWindow(wx.Frame):
atoms = ase.io.read(s, format='xyz') atoms = ase.io.read(s, format='xyz')
cluster_viewer.set_atoms(atoms, rescale=True, center=True) cluster_viewer.set_atoms(atoms, rescale=True, center=True)
cluster_viewer.rotate_atoms(45., 45.) cluster_viewer.rotate_atoms(45., 45.)
cluster_viewer.show_emitter(True) #cluster_viewer.show_emitter(True)
win.Show() win.Show()
def on_viewparameters(self, event): def on_viewparameters(self, event):

20
src/msspec/msspecgui/msspec/gui/clusterviewer.py
View File

@ -93,10 +93,10 @@ class ClusterViewer(wx.Window):
self.Bind(wx.EVT_RIGHT_UP, self.__evt_right_up_cb) self.Bind(wx.EVT_RIGHT_UP, self.__evt_right_up_cb)
self.Bind(wx.EVT_TIMER, self.__evt_timer_cb, self.timer) self.Bind(wx.EVT_TIMER, self.__evt_timer_cb, self.timer)
def show_emitter(self, show=True): def show_emitter(self, show=True, alpha=0.25):
_opts = self.sprites_opts.copy() _opts = self.sprites_opts.copy()
if show: if show:
self.sprites_opts['alpha'] = 0.25 self.sprites_opts['alpha'] = alpha
self.sprites_opts['glow'] = False self.sprites_opts['glow'] = False
else: else:
self.sprites_opts = _opts.copy() self.sprites_opts = _opts.copy()
@ -325,6 +325,19 @@ class ClusterViewer(wx.Window):
self.update_drawing() self.update_drawing()
def __evt_mousewheel_cb(self, event): def __evt_mousewheel_cb(self, event):
if wx.GetKeyState(wx.WXK_CONTROL):
alpha = self.sprites_opts['alpha']
rot = event.GetWheelRotation()
if rot > 0:
alpha *= 1.2
alpha = min(1, alpha)
elif rot < 0:
alpha /= 1.2
alpha = max(0, alpha)
self.sprites_opts['alpha'] = alpha
self.create_atom_sprites()
self.update_drawing()
else:
rot = event.GetWheelRotation() rot = event.GetWheelRotation()
self.timer.Stop() self.timer.Stop()
self.timer.Start(self.refresh_delay) self.timer.Start(self.refresh_delay)
@ -385,7 +398,7 @@ class ClusterViewer(wx.Window):
if glow: if glow:
gradient = cairo.RadialGradient(radius, radius, radius / 2, gradient = cairo.RadialGradient(radius, radius, radius / 2,
radius, radius, radius) radius, radius, radius)
gradient.add_color_stop_rgba(0., 1., 1., 1., .5) gradient.add_color_stop_rgba(0., 1., 1., 1., .5*alpha)
gradient.add_color_stop_rgba(0.5, 1., 1., 1., 0) gradient.add_color_stop_rgba(0.5, 1., 1., 1., 0)
gradient.add_color_stop_rgba(1., 1., 1., 1., 0.) gradient.add_color_stop_rgba(1., 1., 1., 1., 0.)
ctx.set_source(gradient) ctx.set_source(gradient)
@ -463,7 +476,6 @@ class ClusterViewer(wx.Window):
self.__outer_margin *= 1.1 self.__outer_margin *= 1.1
def create_background_sprite(self, w, h): def create_background_sprite(self, w, h):
surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, w, h) surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, w, h)
ctx = cairo.Context(surface) ctx = cairo.Context(surface)

6
src/msspec/parameters.py
View File

@ -536,9 +536,9 @@ class SpecParameters(BaseParameters):
default=3, fmt='d'), default=3, fmt='d'),
Parameter('eigval_ispectrum_ne', types=int, limits=[0, 1], Parameter('eigval_ispectrum_ne', types=int, limits=[0, 1],
default=1, fmt='d'), default=1, fmt='d'),
Parameter('eigval_ipwm', types=int, limits=[-4, 4], default=0, Parameter('eigval_ipwm', types=int, limits=[-4, 4], default=4,
fmt='d'), fmt='d'),
Parameter('eigval_method', types=str, default='AITK', Parameter('eigval_method', types=str, default='EPSI',
allowed_values=['AITK', 'RICH', 'SALZ', 'EPSI', 'EPSG', allowed_values=['AITK', 'RICH', 'SALZ', 'EPSI', 'EPSG',
'RHOA', 'THET', 'LEGE', 'CHEB', 'OVER', 'RHOA', 'THET', 'LEGE', 'CHEB', 'OVER',
'DURB', 'DLEV', 'TLEV', 'ULEV', 'VLEV', 'DURB', 'DLEV', 'TLEV', 'ULEV', 'VLEV',
@ -1026,7 +1026,7 @@ class SourceParameters(BaseParameters):
:ref:`this figure <ped_full_picture>` for questions regarding the proper :ref:`this figure <ped_full_picture>` for questions regarding the proper
orientation. orientation.
""")), """)),
Parameter('phi', types=(int, float), limits=(0., 360.), Parameter('phi', types=(int, float), limits=(-180., 180.),
unit=UREG.degree, default=0., doc=textwrap.dedent(""" unit=UREG.degree, default=0., doc=textwrap.dedent("""
The azimuthal angle of the photon incidence (in degrees). Please refer to The azimuthal angle of the photon incidence (in degrees). Please refer to
:ref:`this figure <ped_full_picture>` for questions regarding the proper :ref:`this figure <ped_full_picture>` for questions regarding the proper

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

@ -8,7 +8,7 @@ DEBUG:=0
includes := -I./memalloc/ -I./cluster_gen/ -I./common_sub -I./renormalization includes := -I./memalloc/ -I./cluster_gen/ -I./common_sub -I./renormalization
includes += -I./phd_se_noso_nosp_nosym includes += -I./phd_se_noso_nosp_nosym
includes += -I./eig/common -I./eig/new_mi -I./eig/new_pw includes += -I./eig/common -I./eig/mi -I./eig/pw
memalloc_src:=memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f memalloc_src:=memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f
memalloc_obj:=$(patsubst %.f,%.o, $(memalloc_src)) memalloc_obj:=$(patsubst %.f,%.o, $(memalloc_src))
@ -25,6 +25,9 @@ renormalization_obj:=$(patsubst %.f,%.o, $(renormalization_src))
phd_se_noso_nosp_nosym_src:=$(filter-out phd_se_noso_nosp_nosym/main.f, $(wildcard phd_se_noso_nosp_nosym/*.f)) phd_se_noso_nosp_nosym_src:=$(filter-out phd_se_noso_nosp_nosym/main.f, $(wildcard phd_se_noso_nosp_nosym/*.f))
phd_se_noso_nosp_nosym_obj:=$(patsubst %.f,%.o, $(phd_se_noso_nosp_nosym_src)) phd_se_noso_nosp_nosym_obj:=$(patsubst %.f,%.o, $(phd_se_noso_nosp_nosym_src))
phd_mi_noso_nosp_nosym_src:=$(filter-out phd_mi_noso_nosp_nosym/main.f, $(wildcard phd_mi_noso_nosp_nosym/*.f))
phd_mi_noso_nosp_nosym_obj:=$(patsubst %.f,%.o, $(phd_mi_noso_nosp_nosym_src))
eig_common_src:=$(wildcard eig/common/*.f) eig_common_src:=$(wildcard eig/common/*.f)
eig_common_obj:=$(patsubst %.f,%.o, $(eig_common_src)) eig_common_obj:=$(patsubst %.f,%.o, $(eig_common_src))
@ -37,9 +40,10 @@ eig_pw_obj:=$(patsubst %.f,%.o, $(eig_pw_src))
objects_src := $(memalloc_src) $(cluster_gen_src) $(common_sub_src) objects_src := $(memalloc_src) $(cluster_gen_src) $(common_sub_src)
objects_src += $(renormalization_src) $(phd_se_noso_nosp_nosym_src) objects_src += $(renormalization_src) $(phd_se_noso_nosp_nosym_src)
objects_src += $(eig_common_src) $(eig_mi_src) $(eig_pw_src) objects_src += $(eig_common_src) $(eig_mi_src) $(eig_pw_src)
objects_src += $(phd_mi_noso_nosp_nosym_src)
objects:=$(patsubst %.f,%.o, $(objects_src)) objects:=$(patsubst %.f,%.o, $(objects_src))
libs_targets := phd_se_noso_nosp_nosym.target eig_mi.target eig_pw.target libs_targets := phd_se_noso_nosp_nosym.target phd_mi_noso_nosp_nosym.target eig_mi.target eig_pw.target
EXE=prog EXE=prog
@ -62,6 +66,11 @@ phd_se_noso_nosp_nosym.target: $(memalloc_obj) $(cluster_gen_obj) $(common_sub_o
@$(F2PY) $(includes) $^ $(F2PY_OPTS) -c -m $(patsubst %.target, %, $@) phd_se_noso_nosp_nosym/main.f @$(F2PY) $(includes) $^ $(F2PY_OPTS) -c -m $(patsubst %.target, %, $@) phd_se_noso_nosp_nosym/main.f
@touch $@ @touch $@
phd_mi_noso_nosp_nosym.target: $(memalloc_obj) $(cluster_gen_obj) $(common_sub_obj) $(renormalization_obj) $(phd_mi_noso_nosp_nosym_obj)
@echo "building Python binding..."
@$(F2PY) $(includes) $^ $(F2PY_OPTS) -c -m $(patsubst %.target, %, $@) phd_mi_noso_nosp_nosym/main.f
@touch $@
eig_mi.target: $(memalloc_obj) $(cluster_gen_obj) $(common_sub_obj) $(renormalization_obj) $(eig_common_obj) $(eig_mi_obj) eig_mi.target: $(memalloc_obj) $(cluster_gen_obj) $(common_sub_obj) $(renormalization_obj) $(eig_common_obj) $(eig_mi_obj)
@echo "building Python binding..." @echo "building Python binding..."
@$(F2PY) $(includes) $^ $(F2PY_OPTS) -c -m $(patsubst %.target, %, $@) eig/mi/main.f @$(F2PY) $(includes) $^ $(F2PY_OPTS) -c -m $(patsubst %.target, %, $@) eig/mi/main.f

196
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/axb_mat_ms_la.f Normal file
View File

@ -0,0 +1,196 @@
C
C=======================================================================
C
SUBROUTINE INV_MAT_MS(JE,TAU)
C
C This subroutine stores the multiple scattering matrix and computes
C the scattering path operator TAU^{j 0} exactly, without explicitely
C using the inverse matrix.
C
C (Photoelectron case)
C
C Last modified : 28 Mar 2007
C
USE DIM_MOD
C
USE COOR_MOD
USE INIT_L_MOD
USE TRANS_MOD
C PARAMETER(NLTWO=2*NL_M)
C
COMPLEX*16 HL1(0:2*NL_M),SM(LINMAX*NATCLU_M,LINMAX*NATCLU_M)
COMPLEX*16 IN(LINMAX*NATCLU_M,LINMAX)
COMPLEX*16 SUM_L,ONEC,IC,ZEROC
COMPLEX*16 YLM(0:2*NL_M,-2*NL_M:2*NL_M),TLJ,TLK,EXPKJ
C
COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M)
C
REAL*8 PI,ATTKJ,GNT(0:N_GAUNT),XKJ,YKJ,ZKJ,RKJ,ZDKJ,KRKJ
C
INTEGER IPIV(LINMAX*NATCLU_M)
C
CHARACTER*1 CH
C
DATA PI /3.1415926535898D0/
C
ONEC=(1.D0,0.D0)
IC=(0.D0,1.D0)
ZEROC=(0.D0,0.D0)
IBESS=3
CH='N'
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)*
1 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
IF(JTYP.EQ.1) THEN
IN(KLIN,JLIN)=ONEC
ENDIF
ELSE
SM(KLIN,JLIN)=-TLK*SUM_L
IF(JTYP.EQ.1) THEN
IN(KLIN,JLIN)=ZEROC
ENDIF
ENDIF
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
LW2=(LMAX(1,JE)+1)*(LMAX(1,JE)+1)
C
C Partial inversion of the multiple scattering matrix MS and
C multiplication by T : the LAPACK subroutine performing
C
C A * x = b
C
C is used where b is the block column corresponding to
C the absorber 0 in the identity matrix. x is then TAU^{j 0}.
C
CALL ZGETRF(JLIN,JLIN,SM,LINMAX*NATCLU_M,IPIV,INFO1)
IF(INFO1.NE.0) THEN
WRITE(6,*) ' ---> INFO1 =',INFO1
ELSE
CALL ZGETRS(CH,JLIN,LW2,SM,LINMAX*NATCLU_M,IPIV,
1 IN,LINMAX*NATCLU_M,INFO)
ENDIF
C
C Storage of the Tau matrix
C
JLIN=0
DO JTYP=1,N_PROT
NBTYPJ=NATYP(JTYP)
LMJ=LMAX(JTYP,JE)
DO JNUM=1,NBTYPJ
JATL=NCORR(JNUM,JTYP)
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)
C
DO LK=0,LMK
ILK=LK*LK+LK+1
DO MK=-LK,LK
INDK=ILK+MK
KLIN=KLIN+1
IF((JATL.EQ.1).AND.(LJ.LE.LF2)) THEN
TAU(INDK,INDJ,KATL)=CMPLX(IN(KLIN,JLIN)*TLJ)
ENDIF
ENDDO
ENDDO
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
ENDDO
ENDDO
C
RETURN
C
END

121
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/coumat.f Normal file
View File

@ -0,0 +1,121 @@
C
C=======================================================================
C
SUBROUTINE COUMAT(ITL,MI,LF,MF,DELTA,RADIAL,MATRIX)
C
C This routine calculates the spin-independent PhD optical matrix
C elements for dipolar excitations. It is stored in
C MATRIX(JDIR,JPOL)
C
C Here, the conventions are :
C
C IPOL=1 : linearly polarized light
C IPOL=2 : circularly polarized light
C
C JPOL=1 : +/x polarization for circular/linear light
C JPOL=2 : -/y polarization for circular/linear light
C
C When IDICHR=0, JDIR = 1,2 and 3 correspond respectively to the x,y
C and z directions for the linear polarization. But for IDICHR=1,
C these basis directions are those of the position of the light.
C
C Last modified : 8 Dec 2008
C
USE DIM_MOD
C
USE INIT_L_MOD , L2 => NNL, L3 => LF1, L4 => LF2, L5 => ISTEP_LF
USE SPIN_MOD , I1 => ISPIN, N1 => NSPIN, N2 => NSPIN2, I2 => ISFLI
&P, I8 => IR_DIA, N3 => NSTEP
USE TYPCAL_MOD , I3 => IPHI, I4 => IE, I5 => ITHETA, I6 => IFTHET,
& I7 => IMOD, I9 => I_CP, I10 => I_EXT
C
COMPLEX MATRIX(3,2),SUM_1,SUM_2,DELTA,YLM(3,-1:1),RADIAL
COMPLEX ONEC,IC,IL,COEF,PROD
C
REAL RLM(1-NL_M:NL_M-1,1-NL_M:NL_M-1,0:NL_M-1),GNT(0:N_GAUNT)
REAL THETA(3),PHI(3)
C
DATA PI4S3,C_LIN,SQR2 /4.188790,1.447202,1.414214/
DATA PIS2 /1.570796/
C
ONEC=(1.,0.)
IC=(0.,1.)
C
IF(INITL.EQ.0) GOTO 2
C
M=MF-MI
C
IF(MOD(LF,4).EQ.0) THEN
IL=ONEC
ELSEIF(MOD(LF,4).EQ.1) THEN
IL=IC
ELSEIF(MOD(LF,4).EQ.2) THEN
IL=-ONEC
ELSEIF(MOD(LF,4).EQ.3) THEN
IL=-IC
ENDIF
C
CALL GAUNT(LI,MI,LF,MF,GNT)
C
IF(ITL.EQ.0) THEN
c COEF=CEXP(IC*DELTA)*CONJG(IL)
COEF=CEXP(IC*DELTA)*IL
ELSE
IF(IDICHR.EQ.0) THEN
c COEF=PI4S3*CONJG(IL)
COEF=PI4S3*IL
ELSE
c COEF=C_LIN*CONJG(IL)
COEF=C_LIN*IL
ENDIF
ENDIF
C
PROD=COEF*RADIAL*GNT(1)
C
IF(IDICHR.EQ.0) THEN
YLM(1,-1)=(0.345494,0.)
YLM(1,0)=(0.,0.)
YLM(1,1)=(-0.345494,0.)
YLM(2,-1)=(0.,-0.345494)
YLM(2,0)=(0.,0.)
YLM(2,1)=(0.,-0.345494)
YLM(3,-1)=(0.,0.)
YLM(3,0)=(0.488602,0.)
YLM(3,1)=(0.,0.)
C
DO JDIR=1,3
MATRIX(JDIR,1)=PROD*CONJG(YLM(JDIR,M))
ENDDO
C
ELSEIF(IDICHR.GE.1) THEN
C
THETA(1)=PIS2
PHI(1)=0.
THETA(2)=PIS2
PHI(2)=PIS2
THETA(3)=0.
PHI(3)=0.
C
DO JDIR=1,3
CALL DJMN(THETA(JDIR),RLM,1)
SUM_1=RLM(-1,M,1)*PROD*CEXP((0.,-1.)*M*PHI(JDIR))
SUM_2=RLM(1,M,1)*PROD*CEXP((0.,-1.)*M*PHI(JDIR))
IF(IPOL.EQ.2) THEN
MATRIX(JDIR,1)=SQR2*SUM_1
MATRIX(JDIR,2)=SQR2*SUM_2
ELSEIF(ABS(IPOL).EQ.1) THEN
MATRIX(JDIR,1)=(SUM_2-SUM_1)
MATRIX(JDIR,2)=(SUM_2+SUM_1)*IC
ENDIF
ENDDO
ENDIF
GOTO 1
C
2 DO JDIR=1,3
MATRIX(JDIR,1)=ONEC
MATRIX(JDIR,2)=ONEC
ENDDO
C
1 RETURN
C
END

85
src/msspec/spec/fortran/phd_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/phd_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/phd_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

5123
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/lapack_axb.f Normal file
View File

File diff suppressed because it is too large Load Diff

20
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/main.f Normal file
View File

@ -0,0 +1,20 @@
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_PHD_NS_MI()
END SUBROUTINE RUN

1683
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/main_phd_ns_mi.f Normal file
View File

File diff suppressed because it is too large Load Diff

1131
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/phddif_mi.f Normal file
View File

File diff suppressed because it is too large Load Diff

106
src/msspec/spec/fortran/phd_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

769
src/msspec/spec/fortran/phd_mi_noso_nosp_nosym/treat_phd.f Normal file
View File

@ -0,0 +1,769 @@
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
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(J
&E),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)*NTHE
&TA*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),S
&R_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),VOLD
&IR_1,VOLDIF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),VO
&LDIR2_1,VOLDIF2_1
ENDIF
WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TOTD
&IR_1,TOTDIF_1
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,3) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TO
&TDIR2_1,TOTDIF2_1
ENDIF
ELSE
WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),VOL
&DIR_1,VOLDIF_1,VOLDIR_2,VOLDIF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JVOL,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),V
&OLDIR2_1,VOLDIF2_1,VOLDIR2_2,VOLDIF2_2
ENDIF
WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),TOT
&DIR_1,TOTDIF_1,TOTDIR_2,TOTDIF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JTOT,DTHETA(JTHETA),DPHI(JPHI2),ECIN(JE),T
&OTDIR2_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),EC
&IN(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(J
&E),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 + JP
&HI
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(J
&E),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),VO
&LDIR_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),TO
&TDIR_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),V
&OLDIR_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),T
&OTDIR_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 + J
&PHI
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),S
&R_1,SF_1,SR_2,SF_2
IF(I_EXT.EQ.-1) THEN
WRITE(IUO2,23) JPLAN,DTHETA(JTHETA),DPHI(JPHI2),ECIN(J
&E),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,E1
&2.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 I
&NCLUDE FILE >>>>>>>>>>',/,8X,'<<<<<<<<<< SHOULD BE AT
&LEAST ',I6,' >>>>>>>>>>')
888 FORMAT(A72)
C
6 RETURN
C
END

335
src/msspec/spec/fortran/phd_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

3
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/do_main.f
View File

@ -564,7 +564,8 @@ C
c CALL PLOTFD_A(A,LMAX_A,ITL_A,NL1_A,NAT2_A,NE_A) c CALL PLOTFD_A(A,LMAX_A,ITL_A,NL1_A,NAT2_A,NE_A)
ENDIF ENDIF
WRITE(IUO1,57) WRITE(IUO1,57)
STOP C STOP
GO TO 999
C C
8 IF(IBAS.EQ.0) THEN 8 IF(IBAS.EQ.0) THEN
C C

4
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths1.f
View File

@ -1,8 +1,8 @@
C C
C======================================================================= C=======================================================================
C C
SUBROUTINE FINDPATHS(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,PHIMI SUBROUTINE FINDPATHS(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
&,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C C
C This routine generates all the paths and filters them according to the C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).

18
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths2.f
View File

@ -1,8 +1,8 @@
C C
C======================================================================= C=======================================================================
C C
SUBROUTINE FINDPATHS2(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,PHIM SUBROUTINE FINDPATHS2(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
&I,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C C
C This routine generates all the paths and filters them according to the C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
@ -24,17 +24,17 @@ C
USE TRANS_MOD USE TRANS_MOD
USE TLDW_MOD USE TLDW_MOD
USE VARIA_MOD USE VARIA_MOD
C
DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M)
DIMENSION JPOS(NDIF_M,3),R(NDIF_M)
C
C
C C
COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK
COMPLEX IC,COMPL1,PW(0:NDIF_M) COMPLEX IC,COMPL1,PW(0:NDIF_M)
COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M) COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M)
COMPLEX YLM1(0:NL_M,-NL_M:NL_M) COMPLEX YLM1(0:NL_M,-NL_M:NL_M)
COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2 COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2
C
DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M)
DIMENSION JPOS(NDIF_M,3),R(NDIF_M)
C
C
C C
DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/ DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/
C C
@ -321,8 +321,8 @@ C
XMAXT=AMAX1(XMAXT,CABS(PW1)) XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO ENDDO
ENDDO ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND) IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
&=0 & IT(ND)=0
ENDIF ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN

8
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths3.f
View File

@ -1,8 +1,8 @@
C C
C======================================================================= C=======================================================================
C C
SUBROUTINE FINDPATHS3(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,PHIM SUBROUTINE FINDPATHS3(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
&I,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C C
C This routine generates all the paths and filters them according to the C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
@ -321,8 +321,8 @@ C
XMAXT=AMAX1(XMAXT,CABS(PW1)) XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO ENDDO
ENDDO ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND) IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
&=0 & IT(ND)=0
ENDIF ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN

8
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths4.f
View File

@ -1,8 +1,8 @@
C C
C======================================================================= C=======================================================================
C C
SUBROUTINE FINDPATHS4(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,PHIM SUBROUTINE FINDPATHS4(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
&I,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C C
C This routine generates all the paths and filters them according to the C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
@ -321,8 +321,8 @@ C
XMAXT=AMAX1(XMAXT,CABS(PW1)) XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO ENDDO
ENDDO ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND) IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
&=0 & IT(ND)=0
ENDIF ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN

8
src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths5.f
View File

@ -1,8 +1,8 @@
C C
C======================================================================= C=======================================================================
C C
SUBROUTINE FINDPATHS5(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,PHIM SUBROUTINE FINDPATHS5(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI,
&I,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU)
C C
C This routine generates all the paths and filters them according to the C This routine generates all the paths and filters them according to the
C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH).
@ -321,8 +321,8 @@ C
XMAXT=AMAX1(XMAXT,CABS(PW1)) XMAXT=AMAX1(XMAXT,CABS(PW1))
ENDDO ENDDO
ENDDO ENDDO
IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND) IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))
&=0 & IT(ND)=0
ENDIF ENDIF
IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32
IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN

12
src/msspec/utils.py
View File

@ -68,9 +68,9 @@ def center_cluster(atoms, invert=False):
atoms.translate(-0.5*cell_vector) atoms.translate(-0.5*cell_vector)
def cut_sphere(atoms, radius): def cut_sphere(atoms, radius, center=(0, 0, 0)):
assert radius >= 0, "Please give a positive radius value" assert radius >= 0, "Please give a positive radius value"
radii = np.linalg.norm(atoms.positions, axis=1) radii = np.linalg.norm(atoms.positions - center, axis=1)
indices = np.where(radii <= radius)[0] indices = np.where(radii <= radius)[0]
return atoms[indices] return atoms[indices]
@ -232,7 +232,8 @@ def cut_plane(atoms, x=None, y=None, z=None):
indices = np.where(list(map(constraint, atoms.positions)))[0] indices = np.where(list(map(constraint, atoms.positions)))[0]
return atoms[indices] return atoms[indices]
def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, planes=0): def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
planes=0, shape='spherical'):
"""Creates and returns a cluster based on an Atoms object and some parameters. """Creates and returns a cluster based on an Atoms object and some parameters.
@ -318,7 +319,12 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, p
assert (radius <= diameter/2), "The number of planes is too high compared to the diameter." assert (radius <= diameter/2), "The number of planes is too high compared to the diameter."
radius = max(radius, diameter/2) radius = max(radius, diameter/2)
if shape in ('spherical'):
cluster = cut_sphere(cluster, radius=radius + eps) # cut a sphere in our cluster with the diameter which is indicate in the parameters cluster = cut_sphere(cluster, radius=radius + eps) # cut a sphere in our cluster with the diameter which is indicate in the parameters
elif shape in ('cylindrical'):
cluster = cut_cylinder(cluster, radius=radius + eps) # cut a sphere in our cluster with the diameter which is indicate in the parameters
else:
raise NameError('Unkknown shape specifier: \"{}\"'.format(shape))
if planes!=0: if planes!=0:
zcut = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4)))[::-1][planes-1] - eps # calculate where to cut to get the right number of planes zcut = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4)))[::-1][planes-1] - eps # calculate where to cut to get the right number of planes