|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +"""Select one of the 3 functions at the end of main()""" |
| 3 | +__author__ = "Konstantin Klementiev, Roman Chernikov" |
| 4 | +__date__ = "21 May 2020" |
| 5 | +import numpy as np |
| 6 | +import matplotlib.pyplot as plt |
| 7 | +import time |
| 8 | + |
| 9 | +# path to xrt: |
| 10 | +import os, sys; sys.path.append(os.path.join('..', '..')) # analysis:ignore |
| 11 | +import xrt.backends.raycing.sources as rs |
| 12 | + |
| 13 | + |
| 14 | +def intensity_in_transverse_plane(energy, theta, psi, I0): |
| 15 | + plt.imshow(I0[len(energy)//2, :, :], |
| 16 | + extent=[theta[0]*1e6, theta[-1]*1e6, psi[0]*1e6, psi[-1]*1e6]) |
| 17 | + |
| 18 | + |
| 19 | +def main(case, icalc, plot): |
| 20 | + if case == 1: |
| 21 | + angleMaxX = 30e-6 # rad |
| 22 | + angleMaxZ = 30e-6 # rad |
| 23 | + energy = np.linspace(3800, 4150, 351) |
| 24 | + theta = np.linspace(-1, 1, 51) * angleMaxX |
| 25 | + psi = np.linspace(-1, 1, 51) * angleMaxZ |
| 26 | + kwargsCommon = dict( |
| 27 | + eE=3.0, eI=0.5, eEpsilonX=0.263, eEpsilonZ=0.008, |
| 28 | + period=18.5, n=108, K=0.52, |
| 29 | + xPrimeMax=theta[-1]*1e3, zPrimeMax=psi[-1]*1e3) |
| 30 | + kwargsXRT = dict( |
| 31 | + betaX=9.539, betaZ=1.982, |
| 32 | + xPrimeMaxAutoReduce=False, zPrimeMaxAutoReduce=False, |
| 33 | +# targetOpenCL='CPU', |
| 34 | + distE='BW') |
| 35 | + elif case == 2: |
| 36 | + angleMaxX = 8.7e-6 # rad |
| 37 | + angleMaxZ = 8.7e-6 # rad |
| 38 | + energy = np.linspace(7250, 7450, 201) |
| 39 | + theta = np.linspace(-1, 1, 801) * angleMaxX * 32 |
| 40 | + psi = np.linspace(-1, 1, 101) * angleMaxZ * 4 |
| 41 | + kwargsCommon = dict( |
| 42 | + eE=1.72, eI=0.3, eEpsilonX=11.371, eEpsilonZ=0.1, |
| 43 | + period=17, n=82, K=1.071) |
| 44 | + kwargsXRT = dict( |
| 45 | + betaX=1.65, betaZ=1, |
| 46 | + xPrimeMax=theta[-1]*1e3, zPrimeMax=psi[-1]*1e3, |
| 47 | + xPrimeMaxAutoReduce=False, zPrimeMaxAutoReduce=False, |
| 48 | + distE='BW') |
| 49 | + elif case == 3: |
| 50 | + angleMaxX = 8.7e-6 # rad |
| 51 | + angleMaxZ = 8.7e-6 # rad |
| 52 | + energy = np.linspace(500, 10000, 1901) |
| 53 | + theta = np.linspace(-1, 1, 801) * angleMaxX * 32 |
| 54 | + psi = np.linspace(-1, 1, 101) * angleMaxZ * 4 |
| 55 | + kwargsCommon = dict( |
| 56 | + eE=1.72, eI=0.3, eEpsilonX=11.371, eEpsilonZ=0.1, |
| 57 | + period=17, n=82, K=1.071) |
| 58 | + kwargsXRT = dict( |
| 59 | + betaX=1.65, betaZ=1, |
| 60 | + xPrimeMax=theta[-1]*1e3, zPrimeMax=psi[-1]*1e3, |
| 61 | + xPrimeMaxAutoReduce=False, zPrimeMaxAutoReduce=False, |
| 62 | + distE='BW') |
| 63 | + |
| 64 | + kwargsURGENT = dict( |
| 65 | + eSigmaX=(kwargsCommon['eEpsilonX']*kwargsXRT['betaX']*1e3)**0.5, |
| 66 | + eSigmaZ=(kwargsCommon['eEpsilonZ']*kwargsXRT['betaZ']*1e3)**0.5, |
| 67 | + eMin=energy[0], eMax=energy[-1], |
| 68 | + eN=len(energy)-1, |
| 69 | + icalc=icalc, |
| 70 | + xPrimeMax=angleMaxX*1e3, zPrimeMax=angleMaxZ*1e3, |
| 71 | + nx=12, nz=12) |
| 72 | + if kwargsURGENT['icalc'] == 3: |
| 73 | + kwargsCommon['eEpsilonX'] = 0. |
| 74 | + kwargsCommon['eEpsilonZ'] = 0. |
| 75 | + |
| 76 | + fig = plt.figure(figsize=(10, 6)) |
| 77 | + ax = fig.add_subplot(111) |
| 78 | + ax.set_xlabel(u'energy (keV)') |
| 79 | + apX = angleMaxX * 2 * 1e6 |
| 80 | + apZ = angleMaxZ * 2 * 1e6 |
| 81 | + ax.set_ylabel(u'flux through {0:.1f}×{1:.1f} µrad² (ph/s/0.1%bw)'.format( |
| 82 | + apX, apZ)) |
| 83 | + axplot = ax.semilogy if 'log-y' in plot else ax.plot |
| 84 | + |
| 85 | + # xrt Undulator |
| 86 | + kwargs = kwargsCommon.copy() |
| 87 | + kwargs.update(kwargsXRT) |
| 88 | + dtheta = theta[1] - theta[0] |
| 89 | + dpsi = psi[1] - psi[0] |
| 90 | + source = rs.Undulator(**kwargs) |
| 91 | + I0x, l1, l2, l3 = source.intensities_on_mesh(energy, theta, psi) |
| 92 | + whereTheta = np.argwhere(abs(theta) <= angleMaxX) |
| 93 | + wherePsi = np.argwhere(abs(psi) <= angleMaxZ) |
| 94 | + fluxX = I0x[:, |
| 95 | + slice(whereTheta.min(), whereTheta.max()+1), |
| 96 | + slice(wherePsi.min(), wherePsi.max()+1)]\ |
| 97 | + .sum(axis=(1, 2)) * dtheta * dpsi |
| 98 | + axplot(energy*1e-3, fluxX, label='xrt') |
| 99 | + # end xrt Undulator |
| 100 | + |
| 101 | + if case == 3: |
| 102 | + eS, fS = np.loadtxt("misc/bessy.dc0.gz", skiprows=2, usecols=(0, 1), |
| 103 | + unpack=True) |
| 104 | + axplot(eS*1e-3, fS, label='Spectra') |
| 105 | + |
| 106 | + # UrgentUndulator |
| 107 | + import xrt.backends.raycing as raycing |
| 108 | + beamLine = raycing.BeamLine(azimuth=0, height=0) |
| 109 | + kwargs = kwargsCommon.copy() |
| 110 | + kwargs.update(kwargsURGENT) |
| 111 | + sourceU = rs.UndulatorUrgent(beamLine, **kwargs) |
| 112 | + I0u, l1, l2, l3 = sourceU.intensities_on_mesh() |
| 113 | + # add the other 3 quadrants, except the x=0 and z=0 lines: |
| 114 | + I0u = np.concatenate((I0u[:, :0:-1, :], I0u), axis=1) |
| 115 | + I0u = np.concatenate((I0u[:, :, :0:-1], I0u), axis=2) |
| 116 | + fluxU = I0u.sum(axis=(1, 2)) * sourceU.dx * sourceU.dz |
| 117 | + urgentEnergy = sourceU.Es |
| 118 | + axplot(urgentEnergy*1e-3, fluxU, label='Urgent') |
| 119 | + # end UrgentUndulator |
| 120 | + |
| 121 | + ax.set_ylim(1e9, 4e13) |
| 122 | + ax.legend() |
| 123 | + fig.savefig(u'flux_case{0}_xrt_UrgentICALC{1}.png'.format(case, icalc)) |
| 124 | + |
| 125 | + # fig = plt.figure() |
| 126 | + # intensity_in_transverse_plane(energy, theta, psi, I0x) |
| 127 | + # fig = plt.figure() |
| 128 | + # intensity_in_transverse_plane(energy, theta, psi, I0u) |
| 129 | + |
| 130 | + |
| 131 | +if __name__ == '__main__': |
| 132 | + t0 = time.time() |
| 133 | + # ICALC=1 non-zero emittance, finite N |
| 134 | + # ICALC=2 non-zero emittance, infinite N |
| 135 | + # ICALC=3 zero emittance, finite N |
| 136 | + |
| 137 | +# main(case=1, icalc=3, plot='lin-y') |
| 138 | +# main(case=1, icalc=2, plot='lin-y') |
| 139 | +# main(case=1, icalc=1, plot='lin-y') |
| 140 | + |
| 141 | +# main(case=2, icalc=3, plot='lin-y') |
| 142 | +# main(case=2, icalc=2, plot='lin-y') |
| 143 | +# main(case=2, icalc=1, plot='lin-y') |
| 144 | + |
| 145 | +# main(case=3, icalc=3, plot='log-y') |
| 146 | +# main(case=3, icalc=2, plot='log-y') |
| 147 | + main(case=3, icalc=1, plot='log-y') |
| 148 | + |
| 149 | + print("Done in {0} s".format(time.time()-t0)) |
| 150 | + plt.show() |
0 commit comments