From ce6088984b37fa65fc30ebd07fcae098806cacdb Mon Sep 17 00:00:00 2001 From: Marc DeGraef Date: Mon, 19 Feb 2024 08:40:49 -0500 Subject: [PATCH] adds new TEMbook folder with EMTBSM and EMSRCBED programs finally catching up on including all the programs from the CTEM text book in one form or another ... --- NamelistTemplates/EMSRCBED.template | 41 + Source/CTEMbook/CMakeLists.txt | 50 + Source/CTEMbook/EMSRCBED.f90 | 69 ++ Source/CTEMbook/EMTBSM.f90 | 333 +++++ Source/EMsoftOOLib/CMakeLists.txt | 1 + Source/EMsoftOOLib/mod_postscript.f90 | 2 +- .../EMsoftOOLib/program_mods/mod_SRCBED.f90 | 1082 +++++++++++++++++ Source/Source.cmake | 1 + Source/TEM/CMakeLists.txt | 1 - resources/templatecodes.txt | 1 + 10 files changed, 1579 insertions(+), 2 deletions(-) create mode 100755 NamelistTemplates/EMSRCBED.template create mode 100644 Source/CTEMbook/CMakeLists.txt create mode 100644 Source/CTEMbook/EMSRCBED.f90 create mode 100644 Source/CTEMbook/EMTBSM.f90 create mode 100644 Source/EMsoftOOLib/program_mods/mod_SRCBED.f90 diff --git a/NamelistTemplates/EMSRCBED.template b/NamelistTemplates/EMSRCBED.template new file mode 100755 index 0000000..2c71c2b --- /dev/null +++ b/NamelistTemplates/EMSRCBED.template @@ -0,0 +1,41 @@ +&SRCBEDlist +! Note: All values in this template file correspond to the program defaults +! +!------------ +! multi-threading parameters +!------------ +! number of threads to be used in multi-threaded mode + nthreads = 6, +!------------ +! microscope parameters +!------------ +! microscope accelerating voltage [kV] + voltage = 200.0, +! camera length [mm] + camlen = 1000.0, +!------------ +! crystal structure and incident beam parameters +!------------ +! string with the filename of the structure data (e.g., Cu.xtal; default = undefined) + xtalname = 'undefined', +! fundamental reciprocal lattice vector of the systematic row + SRG = 1 0 0, +! incident beam direction as a real space direction vector; must be normal to SRG + SRK = 0 0 1, +! range of reflections to be taken into account (-Grange*SRG to +Grange*SRG); number of beams = 2*Grange+1 + Grange = 4, +! up to 10 Laue centers for G at the center of the diffraction disk, defined as k_t/|g| (GLaue = 0.5 means g in Bragg) + numrows = 10, + ktonG = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0, + ! incident beam convergence angle in mrad + convergence = 5.0, +! foil thickness [nm] + thick = 100.0, +!------------ +! various filenames +!------------ +! full filename of the HDF5 output file [relative to EMdatapathname] + outname = 'undefined', +! tiff prefix to be used for each systematic row pattern [relative to EMdatapathname] + tiffprefix = 'undefined', +/ diff --git a/Source/CTEMbook/CMakeLists.txt b/Source/CTEMbook/CMakeLists.txt new file mode 100644 index 0000000..a7bccfc --- /dev/null +++ b/Source/CTEMbook/CMakeLists.txt @@ -0,0 +1,50 @@ +set(APP_DIR "${EMsoftOO_SOURCE_DIR}/Source/CTEMbook") +set(TMPLT_DIR "${EMsoftOO_SOURCE_DIR}/NamelistTemplates") +set(LIB_SEARCH_DIRS ${CMAKE_LIBRARY_OUTPUT_DIRECTORY} ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + +if((EMsoftOO_ENABLE_HDF5_SUPPORT) AND (EMsoftOO_ENABLE_OpenCL_SUPPORT)) + + # The libraries are specifically laid out this way in order to ensure the link line + # has the correct ordering. This seems to be more important on Windows than on macOS/Linux + set(EXE_LINK_LIBRARIES + ${OpenCL_LIBRARY} + clfortran + EMsoftOOLib + EMOpenCLLib + ) + + + # Add_EMsoftOO_Executable(TARGET EMdddSTEM + # SOURCES ${APP_DIR}/EMdddSTEM.f90 + # LINK_LIBRARIES ${EXE_LINK_LIBRARIES} + # TEMPLATE ${TMPLT_DIR}/EMmdSTEM.template + # SOLUTION_FOLDER EMsoftPublic/TEM + # INCLUDE_DIRS ${EMsoftHDFLib_BINARY_DIR} ${EMOpenCLLib_BINARY_DIR} + # ) + +endif() + +if(EMsoftOO_ENABLE_HDF5_SUPPORT) + + GetHDF5LinkLibraries(EMSOFTOO) + set(EXE_LINK_LIBRARIES ${EXE_LINK_LIBRARIES} ${EMSOFTOO_hdf5LinkLibs} ) + + Add_EMsoftOO_Executable(TARGET EMTBSM + SOURCES ${APP_DIR}/EMTBSM.f90 + LINK_LIBRARIES ${EXE_LINK_LIBRARIES} + SOLUTION_FOLDER EMsoftOOPublic/CTEMbook + INSTALL_PROGRAM TRUE + INCLUDE_DIRS ${EMsoftHDFLib_BINARY_DIR} + ) + + Add_EMsoftOO_Executable(TARGET EMSRCBED + SOURCES ${APP_DIR}/EMSRCBED.f90 + LINK_LIBRARIES ${EXE_LINK_LIBRARIES} + TEMPLATE ${TMPLT_DIR}/EMSRCBED.template + SOLUTION_FOLDER EMsoftOOPublic/CTEMbook + INSTALL_PROGRAM TRUE + INCLUDE_DIRS ${EMsoftHDFLib_BINARY_DIR} + ) + + +endif() diff --git a/Source/CTEMbook/EMSRCBED.f90 b/Source/CTEMbook/EMSRCBED.f90 new file mode 100644 index 0000000..c1ca1db --- /dev/null +++ b/Source/CTEMbook/EMSRCBED.f90 @@ -0,0 +1,69 @@ +! ################################################################### +! Copyright (c) 2013-2024, Marc De Graef/Carnegie Mellon University +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without modification, are +! permitted provided that the following conditions are met: +! +! - Redistributions of source code must retain the above copyright notice, this list +! of conditions and the following disclaimer. +! - Redistributions in binary form must reproduce the above copyright notice, this +! list of conditions and the following disclaimer in the documentation and/or +! other materials provided with the distribution. +! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +! of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! ################################################################### + +!-------------------------------------------------------------------------- +! EMsoft: EMSRCBED.f90 +!-------------------------------------------------------------------------- +! +! PROGRAM: EMSRCBED +! +!> @author Marc De Graef, Carnegie Mellon University +! +!> @brief Systematic Row Convergent Beam Electron Diffraction +!> using scattering matrix approach +! +!> @date 06/07/01 MDG 1.0 original (for the CTEMbook) +!> @date 04/18/13 MDG 2.0 rewrite +!> @date 02/16/13 MDG 3.0 complete rewrite for EMsoftOO +!-------------------------------------------------------------------------- +program EMSRCBED + +use mod_kinds +use mod_global +use mod_EMsoft +use mod_SRCBED + +IMPLICIT NONE + + +character(fnlen) :: progname = 'EMSRCBED.f90' +character(fnlen) :: progdesc = 'CTEMbook: Systematic row convergent beam pattern using scattering matrix' + +type(EMsoft_T) :: EMsoft +type(SRCBED_T) :: SRCBED + + +EMsoft = EMsoft_T( progname, progdesc, tpl = (/ 103 /) ) + +! generate a set of systematic row CBED patterns +SRCBED = SRCBED_T( EMsoft%nmldeffile ) + +! perform the computations +call SRCBED%SRCBED( EMsoft, progname ) + +end program EMSRCBED diff --git a/Source/CTEMbook/EMTBSM.f90 b/Source/CTEMbook/EMTBSM.f90 new file mode 100644 index 0000000..a07ba81 --- /dev/null +++ b/Source/CTEMbook/EMTBSM.f90 @@ -0,0 +1,333 @@ +! ################################################################### +! Copyright (c) 2013-2024, Marc De Graef/Carnegie Mellon University +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without modification, are +! permitted provided that the following conditions are met: +! +! - Redistributions of source code must retain the above copyright notice, this list +! of conditions and the following disclaimer. +! - Redistributions in binary form must reproduce the above copyright notice, this +! list of conditions and the following disclaimer in the documentation and/or +! other materials provided with the distribution. +! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +! of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! ################################################################### + +!-------------------------------------------------------------------------- +! EMsoft:EMTBSM.f90 +!-------------------------------------------------------------------------- +! +! PROGRAM: EMTBSM +! +!> @author Marc De Graef, Carnegie Mellon University +! +!> @brief Two-beam scattering matrix program (CTEM book) +! +!> @detail This is a very simplistic implementation to illustrate how the scattering matrix +!> approach can be implemented; this approach is also used in more sophisticated programs. +! +!> @date 4/28/01 MDG 1.0 original +!> @date 5/27/01 MDG 2.0 f90 +!> @date 4/18/13 MDG 3.0 rewrite +!> @date 2/15/24 MDG 4.0 rewrite for EMsoftOO +!-------------------------------------------------------------------------- +program EMTBSM + +use mod_kinds +use mod_global +use mod_EMsoft + +IMPLICIT NONE + +character(fnlen) :: progname = 'EMTBSM.f90' +character(fnlen) :: progdesc = 'CTEMbook: Two-beam computations: comparing scattering matrix and analytical' + +type(EMsoft_T) :: EMsoft + +! initialize +EMsoft = EMsoft_T( progname, progdesc, tpl = (/ 920 /) ) + +! do the computation +call CalcTBSM(EMsoft, progdesc) + +end program + +!-------------------------------------------------------------------------- +! +! SUBROUTINE: CalcTBSM +! +!> @author Marc De Graef, Carnegie Mellon University +! +!> @brief Actual two-beam scattering matrix computation; this does not use the +!> complex data type, but instead explicitly uses real and imaginary parts, +!> just for fun... +! +!> @date 4/28/01 MDG 1.0 original +!> @date 5/27/01 MDG 2.0 f90 +!> @date 4/18/13 MDG 3.0 rewrite +!> @date 2/16/24 MDG 4.0 complete rewrite for EMsoftOO +!-------------------------------------------------------------------------- +subroutine CalcTBSM(EMsoft, progdesc) + +use mod_kinds +use mod_global +use mod_EMsoft +use mod_io +use mod_crystallography +use mod_symmetry +use mod_diffraction +use mod_postscript +use mod_HDFsupport +use mod_memory +use mod_timing +use mod_TB +use mod_misc + +type(EMsoft_T),INTENT(INOUT) :: EMsoft +character(fnlen),INTENT(IN) :: progdesc + +type(Cell_T) :: Cell +type(Spacegroup_T) :: SPG +type(Diffraction_T) :: Diff +type(IO_T) :: Message +type(Memory_T) :: mem +type(Timing_T) :: timer +type(Postscript_T) :: PS +type(gnode) :: rlp + +integer(kind=irg),parameter :: ns=512, nt=512 +integer(kind=irg) :: ind(3),io_int(1), g(3) +real(kind=sgl) :: bg,xg,xgp,xgpz,wmax,tmax, & + Ar(2,2), Ai(2,2), sg, dz, dsg, p(2),q(2), & + BF(512,512,2),DF(512,512,2),It,Is,io_real(1) +real(kind=dbl),allocatable :: phir(:,:),phii(:,:),SMr(:,:,:),SMi(:,:,:) +integer(kind=irg),allocatable :: imaint(:,:) +character(fnlen) :: xtalname + +call openFortranHDFInterface() + +! first get the crystal data and microscope voltage +call Message%ReadValue(' Enter xtal file name : ', xtalname ) +call cell%getCrystalData(xtalname, SPG, EMsoft) +call cell%calcPositions(SPG, 'v') + +call Diff%GetVoltage(cell) +preg = 2.0 * sngl(cRestmass*cCharge/cPlanck**2)*1.0E-18 +call Diff%setrlpmethod('WK') + +! get the reciprocal lattice vector +call Message%printMessage(' Enter diffraction vector :') +call GetIndex(SPG%getSpaceGrouphexset(), g, 'r') + +! some more parameters +call Message%ReadValue(' Enter maximum value of w: ', io_real, 1) +wmax = io_real(1) +call Message%ReadValue(' Maximum foil thickness : ', io_real, 1) +tmax = io_real(1) + +! normal aborption factor +call Diff%CalcUcg(cell, (/ 0, 0, 0 /)) +rlp = Diff%getrlp() +xgpz= rlp%xgp + +! extinction distance and anomalous absorption length +call Diff%CalcUcg(cell, g) +rlp = Diff%getrlp() +xgp = rlp%xgp +xg = rlp%xg +imanum = 0 +dz = tmax/float(nt-1) + +! compute the scattering matrix SM for each orientation +dsg = 2.0*wmax/float(ns-1) +call Message%printMessage('Starting scattering matrix computation') +io_int(1) = 8*ns +call Message%WriteValue(' Number of computations per * = ', io_int, 1, "(I4)") + +! allocate arrays +mem = Memory_T() +call mem%alloc(SMr, (/ 512,2,2 /), 'SMr') +call mem%alloc(SMi, (/ 512,2,2 /), 'SMi') +call mem%alloc(phir, (/ 512,2 /), 'phir') +call mem%alloc(phii, (/ 512,2 /), 'phii') + +timer = Timing_T() +call timer%Time_tick() + + do i=1,ns + sg = (-wmax+dsg*float(i))/xg + call TBCalcSM(Ar,Ai,sg,dz,xg,xgp,xgpz,bg) + SMr(i,1:2,1:2)=Ar(1:2,1:2) + SMi(i,1:2,1:2)=Ai(1:2,1:2) + end do + +! compute intensities for the first row + do i=1,ns + BF(i,1,1)=SMr(i,1,1)**2+SMi(i,1,1)**2 + DF(i,1,1)=SMr(i,2,1)**2+SMi(i,2,1)**2 + +! the first column of SM is the actual initial wavefunction at z=dz + phir(i,1)=SMr(i,1,1) + phir(i,2)=SMr(i,2,1) + phii(i,1)=SMi(i,1,1) + phii(i,2)=SMi(i,2,1) + end do + +! loop to maximum thickness + do l=2,nt + do k=1,ns + do i=1,2 + p(i)=0.D0 + q(i)=0.D0 + do j=1,2 + p(i)=p(i)+SMr(k,i,j)*phir(k,j)-SMi(k,i,j)*phii(k,j) + q(i)=q(i)+SMi(k,i,j)*phir(k,j)+SMr(k,i,j)*phii(k,j) + end do + end do + phir(k,1:2)=p(1:2) + phii(k,1:2)=q(1:2) + end do + if (mod(l,8).eq.0) then + call Message%printMessage('*',"(1A,$)") + end if + do j=1,ns + BF(j,l,1) = phir(j,1)**2+phii(j,1)**2 + DF(j,l,1) = phir(j,2)**2+phii(j,2)**2 + end do + end do + call timer%Time_tock() + call Message%printMessage('done') + io_real(1) = timer%getInterval() + call Message%WriteValue(' Total computation time [s] ', io_real, 1, "(F10.5)") + call Message%printMessage(' Starting direct analytical solution') + +! next, redo the computation, but this time use TBCalcInten + call timer%Time_tick() + do i=1,ns + sg = (-wmax+dsg*float(i))/xg + do j=1,nt + t = float(j)*dz + call TBCalcInten(It,Is,sg,t,xg,xgp,xgpz,bg) + BF(i,j,2) = It + DF(i,j,2) = Is + end do + if (mod(i,8).eq.0) then + call Message%printMessage('*',"(1A,$)") + end if + end do + call timer%Time_tock() + call Message%printMessage('done') + io_real(1) = timer%getInterval() + call Message%WriteValue(' Total computation time [s] ', io_real, 1, "(F10.5)") + +! compare the two computations + bfdiff = 0.0 + dfdiff = 0.0 + do i=1,ns + do j=1,nt + bfdiff = bfdiff + (BF(i,j,1)-BF(i,j,2))**2 + dfdiff = dfdiff + (DF(i,j,1)-DF(i,j,2))**2 + end do + end do + bfdiff = bfdiff/float(ns)/float(nt) + dfdiff = dfdiff/float(ns)/float(nt) + io_real(1) = bfdiff + call Message%WriteValue(' Average difference in BF ', io_real, 1, "(F10.4)") + io_real(1) = dfdiff + call Message%WriteValue(' Average difference in DF ', io_real, 1, "(F10.4)") + call Message%printMessage(' Images computed, preparing for PS output') + +! open PostScript file + call PS%openFile(progdesc, EMsoft) + call PS%setpspage(0) + call PS%newpage(.TRUE.,'Two Beam Rocking Curves') + call PS%setfont(PSfonts(2),0.16) + call PS%text(2.5,8.6,'Scattering Matrix Solution') + call PS%text(2.5,2.1,'Direct Analytical Solution') + call PS%setfont(PSfonts(2),0.10) + call PS%text(0.5,1.8,'Input file') + call PS%text(2.6,1.8,trim(cell%getFileName()) ) + call PS%text(0.5,1.6,'Active reflection') + call PS%textint(2.5,1.6,' ',g(1)) + call PS%textint(2.6,1.6,' ',g(2)) + call PS%textint(2.7,1.6,' ',g(3)) + call PS%text(0.5,1.4,'Extinction distance [nm]') + call PS%textvar(2.5,1.4,' ',xg) + call PS%text(0.5,1.2,'Anomalous absorption length [nm]') + call PS%textvar(2.5,1.2,' ',xgp) + call PS%text(0.5,1.0,'Normal absorption length [nm]') + call PS%textvar(2.5,1.0,' ',xgpz) + call PS%text(0.5,0.8,'Accelerating Voltage [V]') + call PS%textvar(2.5,0.8,' ', sngl(Diff%getV()) ) + call PS%text(0.5,0.6,'Maximum w') + call PS%textvar(2.5,0.6,' ',wmax) + call PS%text(0.5,0.4,'Maximum thickness [nm]') + call PS%textvar(2.5,0.4,' ',tmax) + +! Bright Field image + call mem%alloc(imaint, (/ ns, nt /), 'imaint') + + imanum = 0 + imaint(1:ns,1:nt)=int(255.0*BF(1:ns,1:nt,1)) + x0=0.2 + y0=5.5 + npx=ns + npy=nt + scl=3.0 + call PS%DumpImageDistort(imaint,imanum,x0,y0,npx,npy,scl,scl) + +! Dark Field image + imaint(1:ns,1:nt)=int(255.0*DF(1:ns,1:nt,1)) + x0=3.4 + y0=5.5 + npx=ns + npy=nt + scl=3.0 + call PS%DumpImageDistort(imaint,imanum,x0,y0,npx,npy,scl,scl) + +! Bright Field image + imaint(1:ns,1:nt)=int(255.0*BF(1:ns,1:nt,2)) + x0=0.2 + y0=2.3 + npx=ns + npy=nt + scl=3.0 + call PS%DumpImageDistort(imaint,imanum,x0,y0,npx,npy,scl,scl) + +! Dark Field image + imaint(1:ns,1:nt)=int(255.0*DF(1:ns,1:nt,2)) + x0=3.4 + y0=2.3 + npx=ns + npy=nt + scl=3.0 + call PS%DumpImageDistort(imaint,imanum,x0,y0,npx,npy,scl,scl) + +! close Postscript file + call PS%closeFile() + +! deallocate arrays +call mem%dealloc(SMr, 'SMr') +call mem%dealloc(SMi, 'SMi') +call mem%dealloc(phir, 'phir') +call mem%dealloc(phii, 'phii') +call mem%dealloc(imaint, 'imaint') + +call closeFortranHDFInterface() + +end subroutine CalcTBSM + + diff --git a/Source/EMsoftOOLib/CMakeLists.txt b/Source/EMsoftOOLib/CMakeLists.txt index 63bef79..5bf413d 100644 --- a/Source/EMsoftOOLib/CMakeLists.txt +++ b/Source/EMsoftOOLib/CMakeLists.txt @@ -155,6 +155,7 @@ set(EMsoftOOLib_SRCS ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_HREBSD.f90 ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_CBED.f90 ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_SRdefect.f90 + ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_SRCBED.f90 ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_STEMDCI.f90 ${EMsoftOOLib_Additional_SRCS} diff --git a/Source/EMsoftOOLib/mod_postscript.f90 b/Source/EMsoftOOLib/mod_postscript.f90 index 336de82..4114d2a 100644 --- a/Source/EMsoftOOLib/mod_postscript.f90 +++ b/Source/EMsoftOOLib/mod_postscript.f90 @@ -100,7 +100,7 @@ module mod_postscript ! font-related stuff -character(25),public,parameter :: PSlbl = "Written by MDG, 2001-2020" +character(25),public,parameter :: PSlbl = "Written by MDG, 2001-2024" character(20),public,parameter :: PSfonts(5) = (/"Symbol ", & "Times-Bold ", & "Times-BoldItalic ", & diff --git a/Source/EMsoftOOLib/program_mods/mod_SRCBED.f90 b/Source/EMsoftOOLib/program_mods/mod_SRCBED.f90 new file mode 100644 index 0000000..99709bc --- /dev/null +++ b/Source/EMsoftOOLib/program_mods/mod_SRCBED.f90 @@ -0,0 +1,1082 @@ +! ################################################################### +! Copyright (c) 2013-2024, Marc De Graef Research Group/Carnegie Mellon University +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without modification, are +! permitted provided that the following conditions are met: +! +! - Redistributions of source code must retain the above copyright notice, this list +! of conditions and the following disclaimer. +! - Redistributions in binary form must reproduce the above copyright notice, this +! list of conditions and the following disclaimer in the documentation and/or +! other materials provided with the distribution. +! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +! of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! ################################################################### + +module mod_SRCBED + !! author: MDG + !! version: 1.0 + !! date: 02/16/24 + !! + !! class definition for the EMSRCBED program + +use mod_kinds +use mod_global + +IMPLICIT NONE + +! namelist for the EMSRCBED program +type, public :: SRCBEDNameListType + integer(kind=irg) :: nthreads + integer(kind=irg) :: SRG(3) + integer(kind=irg) :: SRK(3) + integer(kind=irg) :: Grange + integer(kind=irg) :: numrows + real(kind=sgl) :: ktonG(10) + real(kind=sgl) :: camlen + real(kind=sgl) :: voltage + real(kind=sgl) :: convergence + real(kind=sgl) :: thick + character(fnlen) :: outname + character(fnlen) :: xtalname + character(fnlen) :: tiffprefix +end type SRCBEDNameListType + +! class definition +type, public :: SRCBED_T +private + character(fnlen) :: nmldeffile = 'EMSRCBED.nml' + type(SRCBEDNameListType) :: nml + +contains +private + procedure, pass(self) :: readNameList_ + procedure, pass(self) :: getNameList_ + procedure, pass(self) :: SRCBED_ + procedure, pass(self) :: setnthreads_ + procedure, pass(self) :: getnthreads_ + procedure, pass(self) :: setSRG_ + procedure, pass(self) :: getSRG_ + procedure, pass(self) :: setSRK_ + procedure, pass(self) :: getSRK_ + procedure, pass(self) :: setGrange_ + procedure, pass(self) :: getGrange_ + procedure, pass(self) :: setnumrows_ + procedure, pass(self) :: getnumrows_ + procedure, pass(self) :: setktonG_ + procedure, pass(self) :: getktonG_ + procedure, pass(self) :: setvoltage_ + procedure, pass(self) :: getvoltage_ + procedure, pass(self) :: setcamlen_ + procedure, pass(self) :: getcamlen_ + procedure, pass(self) :: setconvergence_ + procedure, pass(self) :: getconvergence_ + procedure, pass(self) :: setthick_ + procedure, pass(self) :: getthick_ + procedure, pass(self) :: setoutname_ + procedure, pass(self) :: getoutname_ + procedure, pass(self) :: setxtalname_ + procedure, pass(self) :: getxtalname_ + procedure, pass(self) :: settiffprefix_ + procedure, pass(self) :: gettiffprefix_ + + generic, public :: getNameList => getNameList_ + generic, public :: readNameList => readNameList_ + generic, public :: SRCBED => SRCBED_ + generic, public :: setnthreads => setnthreads_ + generic, public :: getnthreads => getnthreads_ + generic, public :: setSRG => setSRG_ + generic, public :: getSRG => getSRG_ + generic, public :: setSRK => setSRK_ + generic, public :: getSRK => getSRK_ + generic, public :: setGrange => setGrange_ + generic, public :: getGrange => getGrange_ + generic, public :: setnumrows => setnumrows_ + generic, public :: getnumrows => getnumrows_ + generic, public :: setktonG => setktonG_ + generic, public :: getktonG => getktonG_ + generic, public :: setvoltage => setvoltage_ + generic, public :: getvoltage => getvoltage_ + generic, public :: setcamlen => setcamlen_ + generic, public :: getcamlen => getcamlen_ + generic, public :: setconvergence => setconvergence_ + generic, public :: getconvergence => getconvergence_ + generic, public :: setthick => setthick_ + generic, public :: getthick => getthick_ + generic, public :: setoutname => setoutname_ + generic, public :: getoutname => getoutname_ + generic, public :: setxtalname => setxtalname_ + generic, public :: getxtalname => getxtalname_ + generic, public :: settiffprefix => settiffprefix_ + generic, public :: gettiffprefix => gettiffprefix_ +end type SRCBED_T + +! the constructor routine for this class +interface SRCBED_T + module procedure SRCBED_constructor +end interface SRCBED_T + +contains + +!-------------------------------------------------------------------------- +type(SRCBED_T) function SRCBED_constructor( nmlfile ) result(SRCBED) +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! constructor for the SRCBED_T Class; reads the name list + +IMPLICIT NONE + +character(fnlen), OPTIONAL :: nmlfile + +call SRCBED%readNameList(nmlfile) + +end function SRCBED_constructor + +!-------------------------------------------------------------------------- +subroutine SRCBED_destructor(self) +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! destructor for the SRCBED_T Class + +IMPLICIT NONE + +type(SRCBED_T), INTENT(INOUT) :: self + +call reportDestructor('SRCBED_T') + +end subroutine SRCBED_destructor + +!-------------------------------------------------------------------------- +subroutine readNameList_(self, nmlfile, initonly) +!DEC$ ATTRIBUTES DLLEXPORT :: readNameList_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! read the namelist from an nml file for the SRCBED_T Class + +use mod_io +use mod_EMsoft + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen),INTENT(IN) :: nmlfile + !! full path to namelist file +logical,OPTIONAL,INTENT(IN) :: initonly + !! fill in the default values only; do not read the file + +type(EMsoft_T) :: EMsoft +type(IO_T) :: Message +logical :: skipread = .FALSE. + +integer(kind=irg) :: nthreads +integer(kind=irg) :: SRG(3) +integer(kind=irg) :: SRK(3) +integer(kind=irg) :: Grange +integer(kind=irg) :: numrows +real(kind=sgl) :: ktonG(10) +real(kind=sgl) :: voltage +real(kind=sgl) :: camlen +real(kind=sgl) :: convergence +real(kind=sgl) :: thick +character(fnlen) :: outname +character(fnlen) :: xtalname +character(fnlen) :: tiffprefix + +namelist /SRCBEDlist/ nthreads, SRG, SRK, Grange, numrows, ktonG, voltage, xtalname, convergence, & + thick, outname, tiffprefix, camlen + +nthreads = 6 +voltage = 200.0 +camlen = 1000.0 +SRG = (/ 1, 0, 0 /) +SRK = (/ 0, 0, 1 /) +Grange = 4 +numrows = 10 +ktonG = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /) +convergence = 5.0 +thick = 100.0 +xtalname = 'undefined' +outname = 'undefined' +tiffprefix = 'undefined' + +if (present(initonly)) then + if (initonly) skipread = .TRUE. +end if + +if (.not.skipread) then +! read the namelist file + open(UNIT=dataunit,FILE=trim(nmlfile),DELIM='apostrophe',STATUS='old') + read(UNIT=dataunit,NML=SRCBEDlist) + close(UNIT=dataunit,STATUS='keep') + +! check for required entries + if (trim(xtalname).eq.'undefined') then + call Message%printError('readNameList:',' xtalname file name is undefined in '//nmlfile) + end if + + if (trim(outname).eq.'undefined') then + call Message%printError('readNameList:',' outname file name is undefined in '//nmlfile) + end if + + if (trim(tiffprefix).eq.'undefined') then + call Message%printError('readNameList:',' tiffprefix file name is undefined in '//nmlfile) + end if +end if + +self%nml%nthreads = nthreads +self%nml%SRG = SRG +self%nml%SRK = SRK +self%nml%Grange = Grange +self%nml%numrows = numrows +self%nml%ktonG = ktonG +self%nml%voltage = voltage +self%nml%camlen = camlen +self%nml%convergence = convergence +self%nml%thick = thick +self%nml%outname = outname +self%nml%xtalname = xtalname +self%nml%tiffprefix = tiffprefix + +end subroutine readNameList_ + +!-------------------------------------------------------------------------- +function getNameList_(self) result(nml) +!DEC$ ATTRIBUTES DLLEXPORT :: getNameList_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! pass the namelist for the SRCBED_T Class to the calling program + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +type(SRCBEDNameListType) :: nml + +nml = self%nml + +end function getNameList_ + +!-------------------------------------------------------------------------- +subroutine setnthreads_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setnthreads_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set nthreads in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%nthreads = inp + +end subroutine setnthreads_ + +!-------------------------------------------------------------------------- +function getnthreads_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getnthreads_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get nthreads from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%nthreads + +end function getnthreads_ + +!-------------------------------------------------------------------------- +subroutine setSRG_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setSRG_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set SRG in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp(3) + +self%nml%SRG = inp + +end subroutine setSRG_ + +!-------------------------------------------------------------------------- +function getSRG_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getSRG_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get SRG from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg) :: out(3) + +out = self%nml%SRG + +end function getSRG_ + +!-------------------------------------------------------------------------- +subroutine setSRK_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setSRK_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set SRK in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp(3) + +self%nml%SRK = inp + +end subroutine setSRK_ + +!-------------------------------------------------------------------------- +function getSRK_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getSRK_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get SRK from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg) :: out(3) + +out = self%nml%SRK + +end function getSRK_ + +!-------------------------------------------------------------------------- +subroutine setGrange_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setGrange_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set Grange in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%Grange = inp + +end subroutine setGrange_ + +!-------------------------------------------------------------------------- +function getGrange_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getGrange_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get Grange from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%Grange + +end function getGrange_ + +!-------------------------------------------------------------------------- +subroutine setnumrows_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setnumrows_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set numrows in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%numrows = inp + +end subroutine setnumrows_ + +!-------------------------------------------------------------------------- +function getnumrows_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getnumrows_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get numrows from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%numrows + +end function getnumrows_ + +!-------------------------------------------------------------------------- +subroutine setktonG_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setktonG_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set ktonG in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl), INTENT(IN) :: inp(10) + +self%nml%ktonG = inp + +end subroutine setktonG_ + +!-------------------------------------------------------------------------- +function getktonG_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getktonG_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get ktonG from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl),allocatable :: out(:) + +allocate( out(self%nml%numrows)) +out(1:self%nml%numrows) = self%nml%ktonG(1:self%nml%numrows) + +end function getktonG_ + +!-------------------------------------------------------------------------- +subroutine setvoltage_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setvoltage_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set voltage in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl), INTENT(IN) :: inp + +self%nml%voltage = inp + +end subroutine setvoltage_ + +!-------------------------------------------------------------------------- +function getvoltage_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getvoltage_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get voltage from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl) :: out + +out = self%nml%voltage + +end function getvoltage_ + +!-------------------------------------------------------------------------- +subroutine setcamlen_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setcamlen_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set camlen in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl), INTENT(IN) :: inp + +self%nml%camlen = inp + +end subroutine setcamlen_ + +!-------------------------------------------------------------------------- +function getcamlen_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getcamlen_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get camlen from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl) :: out + +out = self%nml%camlen + +end function getcamlen_ + +!-------------------------------------------------------------------------- +subroutine setconvergence_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setconvergence_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set convergence in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl), INTENT(IN) :: inp + +self%nml%convergence = inp + +end subroutine setconvergence_ + +!-------------------------------------------------------------------------- +function getconvergence_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getconvergence_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get convergence from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl) :: out + +out = self%nml%convergence + +end function getconvergence_ + +!-------------------------------------------------------------------------- +subroutine setthick_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setthick_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set thick in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl), INTENT(IN) :: inp + +self%nml%thick = inp + +end subroutine setthick_ + +!-------------------------------------------------------------------------- +function getthick_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getthick_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get thick from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +real(kind=sgl) :: out + +out = self%nml%thick + +end function getthick_ + +!-------------------------------------------------------------------------- +subroutine setoutname_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setoutname_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set outname in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%outname = trim(inp) + +end subroutine setoutname_ + +!-------------------------------------------------------------------------- +function getoutname_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getoutname_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get outname from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = trim(self%nml%outname) + +end function getoutname_ + +!-------------------------------------------------------------------------- +subroutine setxtalname_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: setxtalname_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set xtalname in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%xtalname = trim(inp) + +end subroutine setxtalname_ + +!-------------------------------------------------------------------------- +function getxtalname_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: getxtalname_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get xtalname from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = trim(self%nml%xtalname) + +end function getxtalname_ + +!-------------------------------------------------------------------------- +subroutine settiffprefix_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: settiffprefix_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! set tiffprefix in the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%tiffprefix = trim(inp) + +end subroutine settiffprefix_ + +!-------------------------------------------------------------------------- +function gettiffprefix_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: gettiffprefix_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! get tiffprefix from the SRCBED_T class + +IMPLICIT NONE + +class(SRCBED_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = trim(self%nml%tiffprefix) + +end function gettiffprefix_ + +!-------------------------------------------------------------------------- +subroutine SRCBED_( self, EMsoft, progname ) +!DEC$ ATTRIBUTES DLLEXPORT :: SRCBED_ +!! author: MDG +!! version: 1.0 +!! date: 02/16/24 +!! +!! perform the computations + +use mod_kinds +use mod_global +use mod_EMsoft +use mod_crystallography +use mod_diffraction +use mod_symmetry +use mod_math +use mod_io +use mod_HDFsupport +use mod_memory +use mod_timing +use mod_initializers +use mod_image + +use, intrinsic :: iso_fortran_env + +IMPLICIT NONE + +class(SRCBED_T),INTENT(INOUT) :: self +type(EMsoft_T),INTENT(INOUT) :: EMsoft +character(fnlen),INTENT(IN) :: progname + +type(Cell_T) :: Cell +type(Spacegroup_T) :: SPG +type(Diffraction_T) :: Diff +type(DynType) :: Dyn +type(IO_T) :: Message +type(Memory_T) :: mem +type(Timing_T) :: timer +type(gnode) :: rlp + +real(kind=dbl) :: laL,z0,alp,thc,thb,omega_c,omega_min,omega_max,dmin, mi, ma, camlen, Upz,& + dom,glen,xgpz, io_real(2),sc,omega,exer,sl,thr,zmax,att,gc,gci, lambda +integer(kind=irg) :: g(3),ira,dpcnt,ppi,io_int(2),nn,izero,npix,i,j,numi,n,l,ll,np2,nps,nsl,& + is,iq,npx,npy,ipos,istart,istop,rowmax,imo,k, wpix,hkl(3), ii, imanum +character(1) :: ans,c +complex(kind=dbl) :: czero = cmplx(0.D0,0.D0,dbl), cone = cmplx(1.D0,0.D0,dbl) +logical :: overlap,np +character(fnlen) :: TIFF_filename +character(3) :: tnum +real(kind=dbl),allocatable :: row(:,:), kt(:) +real(kind=dbl),parameter :: xoff(0:5)=(/0.0,3.3125,0.0,3.3125,0.0,3.3125/),yoff(0:5)=(/6.0,6.0,3.0,3.0,0.0,0.0/) +complex(kind=dbl),allocatable :: SMz(:,:,:),disk(:,:,:),p(:),Sphiz(:,:), DHWMz(:,:) +complex(kind=dbl),allocatable :: q(:,:),qin(:,:),qout(:,:),r(:,:), Az(:,:) + +! declare variables for use in object oriented image module +integer :: iostat +character(len=128) :: iomsg +logical :: isInteger +type(image_t) :: im +integer(int8) :: i8 (3,4) +integer(int8), allocatable :: TIFF_image(:,:) + +call openFortranHDFInterface() + +associate( enl => self%nml ) + +! first get the crystal data and microscope voltage +dmin = 0.05D0 ! should be set in namelist file +call Diff%setV( dble(enl%voltage) ) +call cell%setFileName( enl%xtalname ) +call Initialize_Cell(cell, Diff, SPG, Diff%Dyn, EMsoft, sngl(dmin), noLUT=.TRUE. ) +call cell%calcPositions(SPG, 'v') +call Diff%setrlpmethod('WK') +lambda = Diff%getWaveLength() + +! normal aborption factor +call Diff%CalcUcg(cell, (/ 0, 0, 0 /)) +rlp = Diff%getrlp() +xgpz= rlp%xgp +Upz = rlp%Upmod + +! general parameters +g = enl%SRG +glen = cell%CalcLength(float(g),'r') +z0 = enl%thick +rowmax = enl%numrows +mem = memory_T() +kt = self%getktonG_() + +! various angles in milli-radians +thb = Diff%CalcDiffAngle(cell, g)*0.5D0*1000.D0 +thc = enl%convergence + +io_real(1) = thb +io_real(2) = thc +call Message%WriteValue(' Bragg and convergence angles [mrad] : ', io_real, 2) + +if (thc.ge.thb) call Message%printMessage(' Warning: the diffraction disks appear to overlap...') + +! camera length +camlen = enl%camlen +laL = lambda * camlen +io_real(1) = lambda +call Message%WriteValue(' wavelength [nm] = ', io_real, 1, "(F10.6)") +io_real(1) = camlen +call Message%WriteValue(' L [mm] = ', io_real, 1, "(f10.2)") +io_real(1) = laL +call Message%WriteValue(' camera length lambda*L [mm nm] = ', io_real, 1, "(f10.5)") + +! determine the number of pixels for a diffraction disk +gc = thc/lambda ! radius of disk in nm^-1 + +! scale bar (sc is the conversion factor from nm-1 to inches) +sc = laL/25.4D0 +gci = gc*sc ! disk radius in inches +! define the number of pixels per inch for the computation +ppi = 300 + +! for 600 dpi output we need this many pixels per disk diameter (odd number) +npix = int(ppi * 2.D0 * gci/1000.D0) +if (mod(npix,2).eq.0) npix=npix+1 +io_int(2) = npix + +! next we need the width of the systematic row image; each image row is npix+2 pixels tall +wpix = nint( 2.D0*thb/thc * (npix-1)/2.D0 ) * 2 * enl%Grange + npix +! the last term manages to get the entire outer disks into the image +io_int(1) = wpix +call Message%WriteValue(' Dimensions of single output image : ', io_int, 2, "(I6,' x ',I6)") + +! the number of beams is the same for all patterns +ira = enl%Grange + +! total number of beams +nn = 2*ira+1 +izero = (nn-1)/2+1 +imanum = 0 + +! absorption factor (amplitude) +att = exp(-cPi*z0/xgpz) +io_real(1) = att +call Message%WriteValue(' Absorption factor (amplitude) : ', io_real, 1) + +! allocate dynamical variables +call mem%alloc(SMz, (/ npix,nn,nn /), 'SMz', initval=czero) +call mem%alloc(Sphiz, (/ npix,nn /), 'Sphiz', initval=czero) +call mem%alloc(Az, (/ nn,nn /), 'Az', initval=czero) +call mem%alloc(p, (/ nn /), 'p', initval=czero) +call mem%alloc(DHWMz, (/ nn,nn /), 'DHWMz', initval=czero) +call mem%alloc(q, (/ nn,nn /), 'q', initval=czero) +call mem%alloc(qin, (/ nn,nn /), 'qin', initval=czero) +call mem%alloc(qout, (/ nn,nn /), 'qout', initval=czero) +call mem%alloc(r, (/ nn,nn /), 'r', initval=czero) +call mem%alloc(disk, (/ nn,npix,npix /), 'disk', initval=czero) +call mem%alloc(row, (/ wpix, npix /), 'row', initval = 0.D0) + +allocate(TIFF_image(wpix, npix)) + +! compute the complex DHW matrix +do i=1,nn + do j=1,nn + hkl=(-ira+i-1)*g-(-ira+j-1)*g + if (i.ne.j) then + call Diff%CalcUcg(cell, hkl) + rlp = Diff%getrlp() + DHWMz(i,j) = cPi*cmplx(-aimag(rlp%qg),real(rlp%qg),dbl) + else + DHWMz(i,j) = cmplx(0.0,0.0,dbl) + endif + end do +end do +call Message%printMessage(' Darwin-Howie-Whelan matrix initialized') + +!======================== +! main loop over patterns + do ii = 1, enl%numrows + +! convert k_t to the alp and omega angles + alp = -2.D0*kt(ii)*thb/1000.D0 + omega_c = cPi*0.5D0+alp + omega_min = omega_c - thc/1000.D0 + omega_max = omega_c + thc/1000.D0 + +! step size in omega angle + dom = (omega_max - omega_min)/float(npix-1) + +! loop over all incident beam directions + numi = 0 + do j=1,npix + +! set omega angle + omega = (omega_min+float(j-1)*dom) + do i=1,nn + n = -ira+i-1 +! exer = excitation error + ! exer = -n*glen*cos(omega)-(1.D0-sqrt(1.D0-(n*lambda*glen*sin(omega))**2))/lambda + exer = -n*glen*0.5D0*(2.D0*cos(omega)+n*lambda*glen)/(1.D0+n*lambda*glen*cos(omega)) + DHWMz(i,i)=cmplx(0.D0,2.D0*cPi*exer,dbl) + end do + +! compute the first slice scattering matrix + sl = z0/100.D0 + nsl = 100 + +! make sure slice thickness does not exceed 0.4 nm (somewhat arbitrary) + do while (sl.gt.0.4D0) + sl = sl/2.D0 + nsl = nsl*2 + end do + +! use the MatrixExponential routine from the math module to compute the +! starting scattering matrix. + call MatrixExponential(DHWMz, Az, dble(sl), 'Pade', nn) + + do i=1,nn + SMz(j,i,1:nn) = Az(i,1:nn) + end do + + end do + +! this completes the first slice scattering matrix; next multiply +! it with itself to the desired thickness +! the izero column of SMz is the actual initial wavefunction at z=sl + Sphiz(1:npix,1:nn)=SMz(1:npix,1:nn,izero) + +! loop to maximum thickness + do l=2,nsl + do ll=1,npix + p(1:nn)=cmplx(0.D0,0.D0,dbl) + do i=1,nn + p(i)=p(i)+sum(SMz(ll,i,1:nn)*Sphiz(ll,1:nn)) + end do + Sphiz(ll,1:nn)=p(1:nn) + end do + end do + +! fill the disk variable with amplitudes; fill vertical columns, then apply a mask + do j=1,npix + do i=1,nn + disk(i,j,1:npix) = Sphiz(j,i) + end do + end do + +! next we apply a circular mask (using fourfold symmetry) + np2 = npix/2 + nps = np2**2 + do i=1,npix/2 + is = (i-np2)**2 + do j=1,npix/2 + iq = is+(j-np2)**2 + if (iq.gt.nps) then + do k=1,nn + disk(k,i,j) = czero + disk(k,npix+1-i,j) = czero + disk(k,i,npix+1-j) = czero + disk(k,npix+1-i,npix+1-j) = czero + end do + end if + end do + end do + + disk = disk * cmplx(att,0.D0) + +! the disks may overlap, but we are adding intensities, not amplitudes +! so the overlap does not matter. +! first the zero order beam + row = 0.D0 + ipos = wpix/2+np2+1 + do i=1,npix + row(ipos-i,1:npix) = abs(disk(izero,i,1:npix))**2 + end do + +! negative disks + do j=1,nn/2 + ipos = wpix/2 - np2 - 1 + int(glen*(-ira+j-1)*ppi*sc) + do i=1,npix + if (ipos+i.gt.0) row(ipos+i,1:npix) = row(ipos+i,1:npix) + abs(disk(j,i,1:npix))**2 + end do + end do + +! positive disks + do j=nn/2+2,nn + ipos = wpix/2 - np2 - 1 + int(glen*(-ira+j-1)*ppi*sc) + do i=1,npix + if (ipos+i.lt.wpix) row(ipos+i,1:npix) = row(ipos+i,1:npix) + abs(disk(j,i,1:npix))**2 + end do + end do + +! deallocate the disk variable + where(row.lt.1.D-5) row = 1.D-5 + row = log10(1.D0/row) + rowmax = maxval(row) + row = row/rowmax + row = 1.D0-row + +! write the systematic row image to a tiff file + mi = minval(row) + ma = maxval(row) + + TIFF_image = nint(255.0*(row-mi)/(ma-mi)) + + ! set up the image_t structure + im = image_t(TIFF_image) + if(im%empty()) call Message%printMessage("EMgetADP","failed to convert array to image") + + write (tnum,"(I3.3)") ii + TIFF_filename = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(enl%tiffprefix)//'_'//tnum//'.tiff' + + ! create the file + call im%write(trim(TIFF_filename), iostat, iomsg) ! format automatically detected from extension + if(0.ne.iostat) then + call Message%printMessage(" failed to write image to file : "//iomsg) + else + call Message%printMessage(' systematic row image written to '//trim(TIFF_filename)) + end if + +end do ! main loop over patterns + +! deallocate all arrays + call mem%dealloc(DHWMz, 'DHWMz') + call mem%dealloc(Az, 'Az') + call mem%dealloc(SMz, 'Smz') + call mem%dealloc(Sphiz, 'Sphiz') + call mem%dealloc(p, 'p') + call mem%dealloc(q, 'q') + call mem%dealloc(qin, 'qin') + call mem%dealloc(qout, 'qout') + call mem%dealloc(r, 'r') + call mem%dealloc(row, 'row') + call mem%dealloc(disk, 'disk') + +end associate + +end subroutine SRCBED_ + +end module mod_SRCBED \ No newline at end of file diff --git a/Source/Source.cmake b/Source/Source.cmake index 585d8a4..caca478 100644 --- a/Source/Source.cmake +++ b/Source/Source.cmake @@ -19,6 +19,7 @@ if( ${EMsoftOO_ENABLE_OpenCL_SUPPORT} ) endif() set(MODALITY_DIRS + CTEMbook Demag DictionaryIndexing # GBs diff --git a/Source/TEM/CMakeLists.txt b/Source/TEM/CMakeLists.txt index ba63594..e58819e 100644 --- a/Source/TEM/CMakeLists.txt +++ b/Source/TEM/CMakeLists.txt @@ -92,5 +92,4 @@ if(EMsoftOO_ENABLE_HDF5_SUPPORT) INCLUDE_DIRS ${EMsoftHDFLib_BINARY_DIR} ) - endif() diff --git a/resources/templatecodes.txt b/resources/templatecodes.txt index 91d1f58..c034dd4 100644 --- a/resources/templatecodes.txt +++ b/resources/templatecodes.txt @@ -78,6 +78,7 @@ 100EMoSLERP 101EMgrainviz 102EMCBED2DQC +103EMSRCBED 105EMDIRAM 106EMEBSDpc 107EMEBSDPCA