From 29a133e10d73fc69a72e509f2114fbf1b2184d9d Mon Sep 17 00:00:00 2001 From: Tchoubar Date: Mon, 30 Sep 2024 22:22:02 -0400 Subject: [PATCH] Minor fixes for sim. of exp. --- cpp/src/clients/python/srwlpy.cpp | 4 +- cpp/src/core/sroptel2.cpp | 46 +++++- cpp/src/core/sroptelm.cpp | 5 +- cpp/src/core/srradstr.cpp | 75 ++++++++- cpp/src/core/srradstr.h | 7 +- env/python/srwpy/srwl_bl.py | 238 ++++++++++++++++++++++------- env/python/srwpy/srwl_uti_opt.py | 75 ++++++++- env/python/srwpy/srwlib.py | 25 +++ env/python/srwpy/uti_math_eigen.py | 6 +- 9 files changed, 405 insertions(+), 76 deletions(-) diff --git a/cpp/src/clients/python/srwlpy.cpp b/cpp/src/clients/python/srwlpy.cpp index 4acc5023..bb1a98d0 100644 --- a/cpp/src/clients/python/srwlpy.cpp +++ b/cpp/src/clients/python/srwlpy.cpp @@ -5067,7 +5067,9 @@ static PyObject* srwlpy_ResizeElecFieldMesh(PyObject *self, PyObject *args) //oFunc = PyObject_GetAttrString(oWfr, "allocate"); //END OCTEST - double arPar[] = {0.,1.}; int nPar = 2; double *pPar = arPar; //[0] with or without FFT, [1]==1 means treatment of quadratic term is allowed + //double arPar[] = {0.,1.}; int nPar = 2; double *pPar = arPar; //[0] with or without FFT, [1]==1 means treatment of quadratic term is allowed + //OC22092024 + double arPar[] = { 0.,1.,1.}; int nPar = 3; double *pPar = arPar; //[0] means with (1) or without (0) FFT, [1] means treatment of quadratic term is allowed (1) or not (0), [2] means correction of Re and Im parts of the E-field based on intensity is allowed (1) or not (0) CopyPyListElemsToNumArray(oPar, 'd', pPar, nPar); ProcRes(srwlResizeElecFieldMesh(&wfr, &newMesh, arPar)); diff --git a/cpp/src/core/sroptel2.cpp b/cpp/src/core/sroptel2.cpp index a286a87f..c1b54da6 100644 --- a/cpp/src/core/sroptel2.cpp +++ b/cpp/src/core/sroptel2.cpp @@ -105,6 +105,11 @@ int srTGenOptElem::PropagateRadiationMeth_0(srTSRWRadStructAccessData* pRadAcces pRadDataSingleE->xStep = pRadAccessData->xStep; pRadDataSingleE->zStart = pRadAccessData->zStart; pRadDataSingleE->zStep = pRadAccessData->zStep; + + pRadDataSingleE->xWfrMin = pRadAccessData->xWfrMin; //HG23072024 (fix of a crash at propagation "from waist") + pRadDataSingleE->xWfrMax = pRadAccessData->xWfrMax; + pRadDataSingleE->zWfrMin = pRadAccessData->zWfrMin; + pRadDataSingleE->zWfrMax = pRadAccessData->zWfrMax; } if(pPrevRadDataSingleE != 0) { @@ -121,7 +126,25 @@ int srTGenOptElem::PropagateRadiationMeth_0(srTSRWRadStructAccessData* pRadAcces if(pRadDataSingleE != pRadAccessData) { - if(result = UpdateGenRadStructSliceConstE_Meth_0(pRadDataSingleE, ie, pRadAccessData)) return result; + if(ie > 0) //OC23072024 + { + srTSRWRadStructAccessData &prevRadDataSingleE = vRadSlices[ie - 1]; + if(!gridParamWereModifInSlices) + { + if((pRadDataSingleE->nx != prevRadDataSingleE.nx) || (pRadDataSingleE->xStart != prevRadDataSingleE.xStart) || (pRadDataSingleE->xStep != prevRadDataSingleE.xStep)) gridParamWereModifInSlices = true; + } + if(!gridParamWereModifInSlices) + { + if((pRadDataSingleE->nz != prevRadDataSingleE.nz) || (pRadDataSingleE->zStart != prevRadDataSingleE.zStart) || (pRadDataSingleE->zStep != prevRadDataSingleE.zStep)) gridParamWereModifInSlices = true; + } + } + + //if(result = UpdateGenRadStructSliceConstE_Meth_0(pRadDataSingleE, ie, pRadAccessData)) return result; + //HG23072024 (fix of a crash at propagation "from waist"): + //OC: does "ie < (neOrig - 1)" cast to 0 or 1 only (never 2)? + if(result = UpdateGenRadStructSliceConstE_Meth_0(pRadDataSingleE, ie, pRadAccessData, ie < (neOrig - 1))) return result; //Updating slice data is required even if it is re-interpolated later + + //if(result = UpdateGenRadStructSliceConstE_Meth_0(pRadDataSingleE, ie, pRadAccessData, ie < (neOrig - 1))) return result; //the above doesn't change the transverse grid parameters in *pRadAccessData //vRadSlices.push_back(*pRadDataSingleE); //this automatically calls destructor, which can eventually delete "emulated" structs! @@ -129,9 +152,10 @@ int srTGenOptElem::PropagateRadiationMeth_0(srTSRWRadStructAccessData* pRadAcces srTSRWRadStructAccessData copyRadDataSingleE(*pRadDataSingleE, false); //OC290813 fixing memory leak(?) //this doesn't assume to copy pBaseRadX, pBaseRadZ copyRadDataSingleE.pBaseRadX = copyRadDataSingleE.pBaseRadZ = 0; copyRadDataSingleE.BaseRadWasEmulated = false; vRadSlices.push_back(copyRadDataSingleE); - - if((pRadDataSingleE->nx != pRadAccessData->nx) || (pRadDataSingleE->xStart != pRadAccessData->xStart) || (pRadDataSingleE->xStep != pRadAccessData->xStep)) gridParamWereModifInSlices = true; - if((pRadDataSingleE->nz != pRadAccessData->nz) || (pRadDataSingleE->zStart != pRadAccessData->zStart) || (pRadDataSingleE->zStep != pRadAccessData->zStep)) gridParamWereModifInSlices = true; + + //OC23072024 (commented-out) + //if((pRadDataSingleE->nx != pRadAccessData->nx) || (pRadDataSingleE->xStart != pRadAccessData->xStart) || (pRadDataSingleE->xStep != pRadAccessData->xStep)) gridParamWereModifInSlices = true; + //if((pRadDataSingleE->nz != pRadAccessData->nz) || (pRadDataSingleE->zStart != pRadAccessData->zStart) || (pRadDataSingleE->zStep != pRadAccessData->zStep)) gridParamWereModifInSlices = true; } } @@ -143,6 +167,20 @@ int srTGenOptElem::PropagateRadiationMeth_0(srTSRWRadStructAccessData* pRadAcces {//to test! if(result = ReInterpolateWfrDataOnNewTransvMesh(vRadSlices, pRadDataSingleE, pRadAccessData)) return result; } + else + {//OC23072024 + if((pRadDataSingleE != 0) && (pRadDataSingleE != pRadAccessData)) + { + pRadAccessData->xStart = pRadDataSingleE->xStart; + pRadAccessData->xStep = pRadDataSingleE->xStep; + pRadAccessData->zStart = pRadDataSingleE->zStart; + pRadAccessData->zStep = pRadDataSingleE->zStep; + pRadAccessData->xWfrMin = pRadDataSingleE->xWfrMin; + pRadAccessData->xWfrMax = pRadDataSingleE->xWfrMax; + pRadAccessData->zWfrMin = pRadDataSingleE->zWfrMin; + pRadAccessData->zWfrMax = pRadDataSingleE->zWfrMax; + } + } if((pRadDataSingleE != 0) && (pRadDataSingleE != pRadAccessData)) delete pRadDataSingleE; if((pPrevRadDataSingleE != 0) && (pPrevRadDataSingleE != pRadAccessData)) delete pPrevRadDataSingleE; diff --git a/cpp/src/core/sroptelm.cpp b/cpp/src/core/sroptelm.cpp index e07f554d..404888ec 100644 --- a/cpp/src/core/sroptelm.cpp +++ b/cpp/src/core/sroptelm.cpp @@ -3068,7 +3068,10 @@ int srTGenOptElem::RadResizeCore(srTSRWRadStructAccessData& OldRadAccessData, sr //srwlPrintTime(":RadResizeCore: init variables 1",&start); //if((!RadResizeStruct.DoNotTreatSpherTerm) && WaveFrontTermCanBeTreated(OldRadAccessData)) - if((!RadResizeStruct.doNotTreatSpherTerm()) && WaveFrontTermCanBeTreated(OldRadAccessData)) //OC090311 + //if((!RadResizeStruct.doNotTreatSpherTerm()) && WaveFrontTermCanBeTreated(OldRadAccessData)) //OC090311 + //*************** + if((!RadResizeStruct.doNotTreatSpherTerm()) && WaveFrontTermCanBeTreated(OldRadAccessData) && (OldRadAccessData.Pres < 1)) //OCTEST06262024 (to test resize in ang. repres.) + //*************** { //Added by SY (for profiling?) at parallelizing SRW via OpenMP: //srwlPrintTime(":RadResizeCore: doNotTreatSpherTerm+WaveFrontTermCanBeTreated",&start); diff --git a/cpp/src/core/srradstr.cpp b/cpp/src/core/srradstr.cpp index ea85e065..87e0146c 100644 --- a/cpp/src/core/srradstr.cpp +++ b/cpp/src/core/srradstr.cpp @@ -4655,6 +4655,13 @@ void srTSRWRadStructAccessData::ResizeCoreXZ(SRWLRadMesh& oldMesh, float* pOldRa //bool OrigWfrQuadTermCanBeTreatedAtResizeX = OldRadAccessData.WfrQuadTermCanBeTreatedAtResizeX; //bool OrigWfrQuadTermCanBeTreatedAtResizeZ = OldRadAccessData.WfrQuadTermCanBeTreatedAtResizeZ; + //OC22092024 + bool allowReAndImCorrFromI = true; + if(arPar != 0) + { + if(arPar[2] == 0.) allowReAndImCorrFromI = false; + } + if(allowTreatQuadPhaseTerm && QuadPhaseTermCanBeTreated()) { //NewRadAccessData.WfrQuadTermCanBeTreatedAtResizeX = OldRadAccessData.WfrQuadTermCanBeTreatedAtResizeX; @@ -4834,7 +4841,8 @@ void srTSRWRadStructAccessData::ResizeCoreXZ(SRWLRadMesh& oldMesh, float* pOldRa srTGenOptElem::InterpolFI(InterpolAux02I, xRel, zRel, BufFI, 0); } (*BufFI) *= AuxFI->fNorm; - srTGenOptElem::ImproveReAndIm(BufF, BufFI); + //srTGenOptElem::ImproveReAndIm(BufF, BufFI); + if(allowReAndImCorrFromI) srTGenOptElem::ImproveReAndIm(BufF, BufFI); //OC22092024 //if(FieldShouldBeZeroed) { *BufF = 0.; *(BufF+1) = 0.; } //OC17102021 (commented-out, since this case is treated before these interpolations) *pEX_New = *BufF; *(pEX_New+1) = *(BufF+1); @@ -4844,15 +4852,18 @@ void srTSRWRadStructAccessData::ResizeCoreXZ(SRWLRadMesh& oldMesh, float* pOldRa if(UseLowOrderInterp_PolCompZ) { srTGenOptElem::InterpolF_LowOrder(InterpolAux02, xRel, zRel, BufF, 2); - srTGenOptElem::InterpolFI_LowOrder(InterpolAux02I, xRel, zRel, BufFI, 1); + //srTGenOptElem::InterpolFI_LowOrder(InterpolAux02I, xRel, zRel, BufFI, 1); + if(allowReAndImCorrFromI) srTGenOptElem::InterpolFI_LowOrder(InterpolAux02I, xRel, zRel, BufFI, 1); //OC22092024 } else { srTGenOptElem::InterpolF(InterpolAux02, xRel, zRel, BufF, 2); - srTGenOptElem::InterpolFI(InterpolAux02I, xRel, zRel, BufFI, 1); + //srTGenOptElem::InterpolFI(InterpolAux02I, xRel, zRel, BufFI, 1); + if(allowReAndImCorrFromI) srTGenOptElem::InterpolFI(InterpolAux02I, xRel, zRel, BufFI, 1); //OC22092024 } (*(BufFI+1)) *= (AuxFI+1)->fNorm; - srTGenOptElem::ImproveReAndIm(BufF+2, BufFI+1); + //srTGenOptElem::ImproveReAndIm(BufF+2, BufFI+1); + if(allowReAndImCorrFromI) srTGenOptElem::ImproveReAndIm(BufF+2, BufFI+1); //OC22092024 //if(FieldShouldBeZeroed) { *(BufF+2) = 0.; *(BufF+3) = 0.; } //OC17102021 (commented-out, since this case is treated before these interpolations) *pEZ_New = *(BufF+2); *(pEZ_New+1) = *(BufF+3); @@ -4860,7 +4871,35 @@ void srTSRWRadStructAccessData::ResizeCoreXZ(SRWLRadMesh& oldMesh, float* pOldRa } } } - if(WaveFrontTermWasTreated) TreatQuadPhaseTerm('a', polComp); + //if(WaveFrontTermWasTreated) TreatQuadPhaseTerm('a', polComp); + //OC21092024 + if(WaveFrontTermWasTreated) + {//Treatment of the Quad. Phase Term has to be done on NEW data! + //Setting temporarily all required member variables in this + + //float *pBaseRadX, *pBaseRadZ; + //waveHndl wRad, wRadX, wRadZ; + //int hStateRadX, hStateRadZ; + //double eStep, eStart, xStep, xStart, zStep, zStart; + //long ne, nx, nz; + + pBaseRadX = pNewRadX; + pBaseRadZ = pNewRadZ; + xStep= xStepNew; xStart = xStartNew; + zStep= zStepNew; zStart = zStartNew; + nx = nxNew; + nz = nzNew; + + TreatQuadPhaseTerm('a', polComp); + + //Setting back the member variables to their original values + pBaseRadX = pOldRadX; + pBaseRadZ = pOldRadZ; + xStep = xStepOld; xStart = xStartOld; + zStep = zStepOld; zStart = zStartOld; + nx = oldMesh.nx; + nz = oldMesh.ny; + } } //************************************************************************* @@ -5085,7 +5124,31 @@ void srTSRWRadStructAccessData::ResizeCoreE(SRWLRadMesh& oldMesh, float* pOldRad } } } - if(WaveFrontTermWasTreated) TreatQuadPhaseTerm('r', polComp); + //if(WaveFrontTermWasTreated) TreatQuadPhaseTerm('r', polComp); + //OC22092024 + if(WaveFrontTermWasTreated) + {//Treatment of the Quad. Phase Term has to be done on NEW data! + //Setting temporarily all required member variables in this + + //float *pBaseRadX, *pBaseRadZ; + //waveHndl wRad, wRadX, wRadZ; + //int hStateRadX, hStateRadZ; + //double eStep, eStart, xStep, xStart, zStep, zStart; + //long ne, nx, nz; + + pBaseRadX = pNewRadX; + pBaseRadZ = pNewRadZ; + eStep= eStepNew; eStart = eStartNew; + ne = neNew; + + TreatQuadPhaseTerm('a', polComp); + + //Setting back the member variables to their original values + pBaseRadX = pOldRadX; + pBaseRadZ = pOldRadZ; + eStep = eStepOld; eStart = eStartOld; + ne = oldMesh.ne; + } } //************************************************************************* diff --git a/cpp/src/core/srradstr.h b/cpp/src/core/srradstr.h index 3ef1ac1f..0c69d6d3 100644 --- a/cpp/src/core/srradstr.h +++ b/cpp/src/core/srradstr.h @@ -801,8 +801,11 @@ class srTSRWRadStructAccessData : public CGenObject { char AnglesXAreSmall = (xMagn < CritRatTransvLong*(::fabs(RobsX))); char AnglesZAreSmall = (zMagn < CritRatTransvLong*(::fabs(RobsZ))); - WfrQuadTermCanBeTreatedAtResizeX = (AnglesXAreSmall && RobsXErrIsSmall); - WfrQuadTermCanBeTreatedAtResizeZ = (AnglesZAreSmall && RobsZErrIsSmall); + //WfrQuadTermCanBeTreatedAtResizeX = (AnglesXAreSmall && RobsXErrIsSmall); + //WfrQuadTermCanBeTreatedAtResizeZ = (AnglesZAreSmall && RobsZErrIsSmall); + //OCTEST20092024 (trying to drop he small angle requirement - perhaps it's not really necessary, consider spherical wave e.g.) + WfrQuadTermCanBeTreatedAtResizeX = (RobsXErrIsSmall); + WfrQuadTermCanBeTreatedAtResizeZ = (RobsZErrIsSmall); return (WfrQuadTermCanBeTreatedAtResizeX || WfrQuadTermCanBeTreatedAtResizeZ); } diff --git a/env/python/srwpy/srwl_bl.py b/env/python/srwpy/srwl_bl.py index fed7ba08..3d6c3284 100644 --- a/env/python/srwpy/srwl_bl.py +++ b/env/python/srwpy/srwl_bl.py @@ -155,7 +155,8 @@ def __init__(self, _name='', _e_beam=None, _mag_approx=None, _mag=None, _gsn_bea #------------------------------------------------------------------------ def set_e_beam(self, _e_beam_name='', _e_beam=None, _i=-1, _ens=-1, _emx=-1, _emy=-1, _dr=0, _x=0, _y=0, _xp=0, _yp=0, _e=None, _de=0, _betax=None, _alphax=None, _etax=None, _etaxp=None, _betay=None, _alphay=None, _etay=0, _etayp=0, - _sigx=None, _sigxp=None, _mxxp=None, _sigy=None, _sigyp=None, _myyp=None): + _sigx=None, _sigxp=None, _mxxp=None, _sigy=None, _sigyp=None, _myyp=None, _z=None): #OC23092024 (allowed setting z-coordinate of initial conditions for the case of presence of quads) + #_sigx=None, _sigxp=None, _mxxp=None, _sigy=None, _sigyp=None, _myyp=None): """Setup Electron Beam :param _e_beam_name: e-beam unique name, e.g. 'NSLS-II Low Beta Day 1' (see srwl_uti_src.py) @@ -292,6 +293,7 @@ def check_positive(d): if(_y is not None): self.eBeam.partStatMom1.y = _y if(_xp is not None): self.eBeam.partStatMom1.xp = _xp if(_yp is not None): self.eBeam.partStatMom1.yp = _yp + if(_z is not None): self.eBeam.partStatMom1.z = _z #OC23092024 if((_dr is not None) and (_dr != 0)): self.eBeam.drift(_dr) #OC16032017 (see above) @@ -878,12 +880,19 @@ def calc_mag_fld(self, _z, _rz, _nz, _x=0, _rx=0, _nx=1, _y=0, _ry=0, _ny=1, _ma if(self.mag_approx is None): raise Exception('Approximate Magnetic Field is not defined') elif(_mag_type == 2): if(self.mag is None): raise Exception('Magnetic Field is not defined') + elif(_mag_type == 3): #OC07092024 + if((self.mag_approx is None) and (self.mag is None)): raise Exception('Magnetic Field is not defined') else: raise Exception('Incorrect Magnetic Field type identificator') magToUse = self.mag_approx if(_mag_type == 2): magToUse = self.mag #print('Using tabulated magnetic field...') + elif(_mag_type == 3): #OC07092024 + if magToUse is None: magToUse = self.mag + elif self.mag is not None: + magToUse = deepcopy(self.mag_approx) + magToUse.combine(self.mag) nTot = _nx*_ny*_nz aux = [0]*nTot @@ -915,6 +924,7 @@ def calc_el_trj(self, _ctst, _ctfi, _np=50000, _mag_type=1, _fname=''): :param _mag_type: "type" of magnetic field to use: 1- "Approximate", referenced by self.mag_approx; 2- "Accurate" (tabulated), referenced by self.mag; + 3- "Combined", i.e. referenced by self.mag and self.mag_approx (test); :param _fname: name of file to save the resulting data to (for the moment, in ASCII format) :return: trajectory structure """ @@ -925,12 +935,19 @@ def calc_el_trj(self, _ctst, _ctfi, _np=50000, _mag_type=1, _fname=''): if(self.mag_approx is None): raise Exception('Approximate Magnetic Field is not defined') elif(_mag_type == 2): if(self.mag is None): raise Exception('Magnetic Field is not defined') + elif(_mag_type == 3): #OC07092024 + if((self.mag_approx is None) and (self.mag is None)): raise Exception('Magnetic Field is not defined') else: raise Exception('Incorrect Magnetic Field type identificator') magToUse = self.mag_approx if(_mag_type == 2): magToUse = self.mag #print('Using tabulated magnetic field...') + elif(_mag_type == 3): #OC07092024 + if magToUse is None: magToUse = self.mag + elif self.mag is not None: + magToUse = deepcopy(self.mag_approx) + magToUse.combine(self.mag) partTraj = SRWLPrtTrj() partTraj.partInitCond = self.eBeam.partStatMom1 @@ -1042,6 +1059,7 @@ def calc_sr_se(self, _mesh, _samp_fact=-1, _meth=2, _rel_prec=0.01, _pol=6, _int :param _mag_type: "type" of magnetic field to use: 1- "Approximate", referenced by self.mag_approx; 2- "Accurate" (tabulated), referenced by self.mag; + 3- "Combined", i.e. referenced by self.mag and self.mag_approx (test); :param _fname: name of file to save the resulting data to (for the moment, in ASCII format) :param _det: detector (instance of SRWLDet) :param _zi: initial lonngitudinal position [m] of electron trajectory for SR calculation @@ -1073,12 +1091,19 @@ def calc_sr_se(self, _mesh, _samp_fact=-1, _meth=2, _rel_prec=0.01, _pol=6, _int if(self.mag_approx is None): raise Exception('Approximate Magnetic Field is not defined') elif(_mag_type == 2): if(self.mag is None): raise Exception('Magnetic Field is not defined') + elif(_mag_type == 3): #OC07092024 + if((self.mag_approx is None) and (self.mag is None)): raise Exception('Magnetic Field is not defined') else: raise Exception('Incorrect Magnetic Field type identificator') magToUse = self.mag_approx if(_mag_type == 2): magToUse = self.mag #print('Using tabulated magnetic field...') + elif(_mag_type == 3): #OC07092024 + if magToUse is None: magToUse = self.mag + elif self.mag is not None: + magToUse = deepcopy(self.mag_approx) + magToUse.combine(self.mag) wfr = SRWLWfr() wfr.allocate(_mesh.ne, _mesh.nx, _mesh.ny) #Numbers of points vs Photon Energy, Horizontal and Vertical Positions @@ -1638,10 +1663,18 @@ def calc_pow_den(self, _mesh, _prec=1, _meth=1, _z_start=0., _z_fin=0., _mag_typ if(self.mag_approx is None): raise Exception('Approximate Magnetic Field is not defined') elif(_mag_type == 2): if(self.mag is None): raise Exception('Magnetic Field is not defined') + elif(_mag_type == 3): #OC07092024 + if((self.mag_approx is None) and (self.mag is None)): raise Exception('Magnetic Field is not defined') else: raise Exception('Incorrect Magnetic Field type identificator') magToUse = self.mag_approx - if(_mag_type == 2): magToUse = self.mag + if(_mag_type == 2): + magToUse = self.mag + elif(_mag_type == 3): #OC07092024 + if magToUse is None: magToUse = self.mag + elif self.mag is not None: + magToUse = deepcopy(self.mag_approx) + magToUse.combine(self.mag) arPrecPar = [_prec, _meth, _z_start, _z_fin, 20000] @@ -1941,7 +1974,7 @@ def calc_wfr_prop(self, _wfr, _pres_ang=0, _pol=6, _int_type=0, _dep_type=3, _fn if(isinstance(_wfr, SRWLWfr) == False): raise Exception('Incorrect wavefront (SRWLWfr) structure') - from datetime import datetime #OCTEST + #from datetime import datetime #OCTEST print('Propagation ... ', end='') t0 = time.time(); @@ -2116,7 +2149,13 @@ def calc_wfr_emit_prop_me(self, _mesh, _sr_samp_fact=1, _sr_meth=2, _sr_rel_prec #print(self.mag) #DEBUG magToUse = self.mag_approx - if(_mag_type == 2): magToUse = self.mag + if(_mag_type == 2): + magToUse = self.mag + elif(_mag_type == 3): #OC08092024 + if magToUse is None: magToUse = self.mag + elif self.mag is not None: + magToUse = deepcopy(self.mag_approx) + magToUse.combine(self.mag) #if((magToUse is None) and (self.gsnBeam is not None)): magToUse = self.gsnBeam #OC15092017 (because _mag is used for SRWLGsnBm when doing partially-coherent simulations in the scope of Gaussian-Schell model) if(magToUse is None): @@ -2801,7 +2840,8 @@ def calc_all(self, _v, _op=None): #16122018 _mxxp=_v.ebm_mxxp, _sigy=_v.ebm_sigy, _sigyp=_v.ebm_sigyp, - _myyp=_v.ebm_myyp) + _myyp=_v.ebm_myyp, #) + _z=_v.ebm_z) #OC23092024 #print('e-beam was set up') @@ -2997,28 +3037,113 @@ def calc_all(self, _v, _op=None): #16122018 #_v.w_mag = 2 #_v.tr_mag = 2 - #---setup magnetic field of a dipole magnet - if hasattr(_v, 'mag_bx') == False: _v.mag_bx = 0 - if hasattr(_v, 'mag_by') == False: _v.mag_by = 0 - if hasattr(_v, 'mag_gn') == False: _v.mag_gn = 0 - if hasattr(_v, 'mag_gs') == False: _v.mag_gs = 0 - if hasattr(_v, 'mag_len') == False: _v.mag_len = 1.5 #? - if hasattr(_v, 'mag_led') == False: _v.mag_led = 0 - if hasattr(_v, 'mag_r') == False: _v.mag_r = 0 - if hasattr(_v, 'mag_zc') == False: _v.mag_zc = 0 - if((_v.mag_bx != 0) or (_v.mag_by != 0) or (_v.mag_gn != 0) or (_v.mag_gs != 0)): - print('Before setting up multipole') #DEBUG - self.set_mag_multipole(#setup dipole / quad magnet parameters - _bx = _v.mag_bx, - _by = _v.mag_by, - _gn = _v.mag_gn, - _gs = _v.mag_gs, - _len = _v.mag_len, - _led = _v.mag_led, - _r = _v.mag_r, - _zc = _v.mag_zc) - #self.mag = None #OC16122016 (commented-out) - + #---setup magnetic field of a dipole/multipole magnet + if(hasattr(_v, 'mag_bx') or hasattr(_v, 'mag_by') or hasattr(_v, 'mag_gn') or hasattr(_v, 'mag_gs')): #OC07092024 + if hasattr(_v, 'mag_bx') == False: _v.mag_bx = 0 + if hasattr(_v, 'mag_by') == False: _v.mag_by = 0 + if hasattr(_v, 'mag_gn') == False: _v.mag_gn = 0 + if hasattr(_v, 'mag_gs') == False: _v.mag_gs = 0 + if hasattr(_v, 'mag_len') == False: _v.mag_len = 1. #? + if hasattr(_v, 'mag_led') == False: _v.mag_led = 0 + if hasattr(_v, 'mag_r') == False: _v.mag_r = 0 + if hasattr(_v, 'mag_zc') == False: _v.mag_zc = 0 + if((_v.mag_bx != 0) or (_v.mag_by != 0) or (_v.mag_gn != 0) or (_v.mag_gs != 0)): + #print('Before setting up multipole') #DEBUG + self.set_mag_multipole(#setup dipole / quad magnet parameters + _bx = _v.mag_bx, + _by = _v.mag_by, + _gn = _v.mag_gn, + _gs = _v.mag_gs, + _len = _v.mag_len, + _led = _v.mag_led, + _r = _v.mag_r, + _zc = _v.mag_zc) + #self.mag = None #OC16122016 (commented-out) + + if(hasattr(_v, 'mag2_bx') or hasattr(_v, 'mag2_by') or hasattr(_v, 'mag2_gn') or hasattr(_v, 'mag2_gs')): #OC07092024 + if not hasattr(_v, 'mag2_bx'): _v.mag2_bx = 0 + if not hasattr(_v, 'mag2_by'): _v.mag2_by = 0 + if not hasattr(_v, 'mag2_gn'): _v.mag2_gn = 0 + if not hasattr(_v, 'mag2_gs'): _v.mag2_gs = 0 + if not hasattr(_v, 'mag2_len'): _v.mag2_len = 1. #? + if not hasattr(_v, 'mag2_led'): _v.mag2_led = 0 + if not hasattr(_v, 'mag2_r'): _v.mag2_r = 0 + if not hasattr(_v, 'mag2_zc'): _v.mag2_zc = 0 + if((_v.mag2_bx != 0) or (_v.mag2_by != 0) or (_v.mag2_gn != 0) or (_v.mag2_gs != 0)): + self.set_mag_multipole(#setup dipole / quad magnet parameters + _bx = _v.mag2_bx, + _by = _v.mag2_by, + _gn = _v.mag2_gn, + _gs = _v.mag2_gs, + _len = _v.mag2_len, + _led = _v.mag2_led, + _r = _v.mag2_r, + _zc = _v.mag2_zc, + _add = 1) + + if(hasattr(_v, 'mag3_bx') or hasattr(_v, 'mag3_by') or hasattr(_v, 'mag3_gn') or hasattr(_v, 'mag3_gs')): #OC07092024 + if not hasattr(_v, 'mag3_bx'): _v.mag3_bx = 0 + if not hasattr(_v, 'mag3_by'): _v.mag3_by = 0 + if not hasattr(_v, 'mag3_gn'): _v.mag3_gn = 0 + if not hasattr(_v, 'mag3_gs'): _v.mag3_gs = 0 + if not hasattr(_v, 'mag3_len'): _v.mag3_len = 1. #? + if not hasattr(_v, 'mag3_led'): _v.mag3_led = 0 + if not hasattr(_v, 'mag3_r'): _v.mag3_r = 0 + if not hasattr(_v, 'mag3_zc'): _v.mag3_zc = 0 + if((_v.mag3_bx != 0) or (_v.mag3_by != 0) or (_v.mag3_gn != 0) or (_v.mag3_gs != 0)): + self.set_mag_multipole(#setup dipole / quad magnet parameters + _bx = _v.mag3_bx, + _by = _v.mag3_by, + _gn = _v.mag3_gn, + _gs = _v.mag3_gs, + _len = _v.mag3_len, + _led = _v.mag3_led, + _r = _v.mag3_r, + _zc = _v.mag3_zc, + _add = 1) + + if(hasattr(_v, 'mag4_bx') or hasattr(_v, 'mag4_by') or hasattr(_v, 'mag4_gn') or hasattr(_v, 'mag4_gs')): #OC07092024 + if not hasattr(_v, 'mag4_bx'): _v.mag4_bx = 0 + if not hasattr(_v, 'mag4_by'): _v.mag4_by = 0 + if not hasattr(_v, 'mag4_gn'): _v.mag4_gn = 0 + if not hasattr(_v, 'mag4_gs'): _v.mag4_gs = 0 + if not hasattr(_v, 'mag4_len'): _v.mag4_len = 1. #? + if not hasattr(_v, 'mag4_led'): _v.mag4_led = 0 + if not hasattr(_v, 'mag4_r'): _v.mag4_r = 0 + if not hasattr(_v, 'mag4_zc'): _v.mag4_zc = 0 + if((_v.mag4_bx != 0) or (_v.mag4_by != 0) or (_v.mag4_gn != 0) or (_v.mag4_gs != 0)): + self.set_mag_multipole(#setup dipole / quad magnet parameters + _bx = _v.mag4_bx, + _by = _v.mag4_by, + _gn = _v.mag4_gn, + _gs = _v.mag4_gs, + _len = _v.mag4_len, + _led = _v.mag4_led, + _r = _v.mag4_r, + _zc = _v.mag4_zc, + _add = 1) + + if(hasattr(_v, 'mag5_bx') or hasattr(_v, 'mag5_by') or hasattr(_v, 'mag5_gn') or hasattr(_v, 'mag5_gs')): #OC07092024 + if not hasattr(_v, 'mag5_bx'): _v.mag5_bx = 0 + if not hasattr(_v, 'mag5_by'): _v.mag5_by = 0 + if not hasattr(_v, 'mag5_gn'): _v.mag5_gn = 0 + if not hasattr(_v, 'mag5_gs'): _v.mag5_gs = 0 + if not hasattr(_v, 'mag5_len'): _v.mag5_len = 1. #? + if not hasattr(_v, 'mag5_led'): _v.mag5_led = 0 + if not hasattr(_v, 'mag5_r'): _v.mag5_r = 0 + if not hasattr(_v, 'mag5_zc'): _v.mag5_zc = 0 + if((_v.mag5_bx != 0) or (_v.mag5_by != 0) or (_v.mag5_gn != 0) or (_v.mag5_gs != 0)): + self.set_mag_multipole(#setup dipole / quad magnet parameters + _bx = _v.mag5_bx, + _by = _v.mag5_by, + _gn = _v.mag5_gn, + _gs = _v.mag5_gs, + _len = _v.mag5_len, + _led = _v.mag5_led, + _r = _v.mag5_r, + _zc = _v.mag5_zc, + _add = 1) + #---setup arbitrary magnetic field source from one file of tabulated 3D magnetic field data if hasattr(_v, 'mag_fn'): if(len(_v.mag_fn) > 0): @@ -3800,38 +3925,40 @@ def calc_all(self, _v, _op=None): #16122018 ['Horizontal Position', 'Vertical Position', sValLabel + ' After Elem. #' + repr(curIntData[0])], ['m', 'm', sValUnit], True) + + plotOK = True #OC12092024 - #if _v.wm and (len(_v.wm_pl) > 0): #OC19012021 - # if (_v.wm_pl == 'xy') or (_v.wm_pl == 'yx') or (_v.wm_pl == 'XY') or (_v.wm_pl == 'YX'): - # #print('2D plot panel is to be prepared') + if _v.wm and (len(_v.wm_pl) > 0): #OC19012021 + if (_v.wm_pl == 'xy') or (_v.wm_pl == 'yx') or (_v.wm_pl == 'XY') or (_v.wm_pl == 'YX'): + #print('2D plot panel is to be prepared') - # sValLabel = 'Flux per Unit Surface' - # sValUnit = 'ph/s/.1%bw/mm^2' - # if(_v.w_u == 0): - # sValLabel = 'Intensity' - # sValUnit = 'a.u.' - # elif(_v.w_u == 2): - # if(_v.w_ft == 't'): - # sValLabel = 'Power Density' - # sValUnit = 'W/mm^2' - # elif(_v.w_ft == 'f'): - # sValLabel = 'Spectral Fluence' - # sValUnit = 'J/eV/mm^2' - - # uti_plot2d1d( - # #int_w0, - # _v.wm_res.arS, - # [_v.wm_res.mesh.xStart, _v.wm_res.mesh.xFin, _v.wm_res.mesh.nx], - # [_v.wm_res.mesh.yStart, _v.wm_res.mesh.yFin, _v.wm_res.mesh.ny], - # 0, #0.5*(mesh_si.xStart + mesh_si.xFin), - # 0, #0.5*(mesh_si.yStart + mesh_si.yFin), - # ['Horizontal Position', 'Vertical Position', sValLabel + ' After Propagation'], - # ['m', 'm', sValUnit], - # True) + sValLabel = 'Flux per Unit Surface' + sValUnit = 'ph/s/.1%bw/mm^2' + if(_v.w_u == 0): + sValLabel = 'Intensity' + sValUnit = 'a.u.' + elif(_v.w_u == 2): + if(_v.w_ft == 't'): + sValLabel = 'Power Density' + sValUnit = 'W/mm^2' + elif(_v.w_ft == 'f'): + sValLabel = 'Spectral Fluence' + sValUnit = 'J/eV/mm^2' + + uti_plot2d1d( + #int_w0, + _v.wm_res.arS, + [_v.wm_res.mesh.xStart, _v.wm_res.mesh.xFin, _v.wm_res.mesh.nx], + [_v.wm_res.mesh.yStart, _v.wm_res.mesh.yFin, _v.wm_res.mesh.ny], + 0, #0.5*(mesh_si.xStart + mesh_si.xFin), + 0, #0.5*(mesh_si.yStart + mesh_si.yFin), + ['Horizontal Position', 'Vertical Position', sValLabel + ' After Propagation'], + ['m', 'm', sValUnit], + True) + plotOK = True #to continue here - plotOK = True if(plotOK): uti_plot_show() @@ -4199,7 +4326,8 @@ def srwl_uti_std_options(): ['wm_fncm', 's', '', 'file name of input coherent modes; if this file name is supplied, the eventual partially-coherent radiation propagation simulation will be done based on propagation of the coherent modes from that file.'], #OC02072021 ['wm_fbk', '', '', 'create backup file(s) with propagated multi-e intensity distribution vs horizontal and vertical position and other radiation characteristics', 'store_true'], - ['wm_pl', 's', 'xy', 'plot the propagated radiaiton intensity distributions in graph(s): ""- dont plot, "x"- vs horizontal position, "y"- vs vertical position, "xy"- vs horizontal and vertical position'], + ['wm_pl', 's', '', 'plot the propagated radiaiton intensity distributions in graph(s): ""- dont plot, "x"- vs horizontal position, "y"- vs vertical position, "xy"- vs horizontal and vertical position'], #OC25072024 + #['wm_pl', 's', 'xy', 'plot the propagated radiaiton intensity distributions in graph(s): ""- dont plot, "x"- vs horizontal position, "y"- vs vertical position, "xy"- vs horizontal and vertical position'], #['ws_fn', 's', '', 'file name for saving single-e (/ fully coherent) wavefront data'], #['wm_fn', 's', '', 'file name for saving multi-e (/ partially coherent) wavefront data'], diff --git a/env/python/srwpy/srwl_uti_opt.py b/env/python/srwpy/srwl_uti_opt.py index 5380abee..389e511f 100644 --- a/env/python/srwpy/srwl_uti_opt.py +++ b/env/python/srwpy/srwl_uti_opt.py @@ -8,12 +8,19 @@ from numpy import fft from array import * from scipy.interpolate import interp2d +from array import * #OC15052024 -try: - from srwlib import * -except: - from oasys_srw.srwlib import * # get oasys_srw: https://github.com/oasys-kit/OASYS1-srwpy +#try: +# from srwlib import * +#except: +# from oasys_srw.srwlib import * # get oasys_srw: https://github.com/oasys-kit/OASYS1-srwpy +try: #OC15052024 + from .srwlib import * + from . import uti_math +except: + from srwlib import * + import uti_math ###################################################################### ## Optical Elements @@ -400,7 +407,6 @@ def srw_uti_mtrl_Param_psd2D(sigma, exponent, PixelWidth, m , n, qr=0, C=None): ## Surface generation from PSD ###################################################################### - def srw_uti_mtrl_psd2D_Prof(Cq, qx, qy, symmetry=False, dist=0, seed=None): ''' parameters (in SI units) @@ -585,4 +591,61 @@ def srw_uti_mtrl_Param_Prof(sigma, exponent, PixelWidth, m , n, qr=0, symmetry=F If the axis direction is defined to be the z axis, and the meridional plane is the y-z plane, sagittal rays intersect the pupil at yp=0. The principal ray is both sagittal and meridional. All other sagittal rays are skew rays. -""" \ No newline at end of file +""" + +###################################################################### +## Misc. Auxiliary function to calculate spectral transmission of beamlines with grating mono from avg. monochromatic intensity distribution and dispersion rate at an exit slit +###################################################################### +def srwl_uti_opt_spec_transm(_ar, _x_grid, _y_grid, _dx, _dy, _re, _ne, _disp_rate, _disp_pl='v', _e0=0., _x0=0., _y0=0., _nx_int=0, _ny_int=0): #OC15052024 + """Auxiliary function to calculate spectral transmission of beamlines with grating mono from avg. monochromatic intensity distribution and dispersion rate at an exit slit + + Args: + _ar_int (python array): input 2D flat array of the monochromatic intensity distribution at the exit slit [ph/s/.1%bw/mm^2] + _x_grid (list/array): list/array specifying grid vs one dimensions (_x_grid[0] is start, _x_grid[1] is end, _x_grid[2] is number of points) + _y_grid (list/array): list/array specifying grid vs another dimensions (_y_grid[0] is start, _y_grid[1] is end, _y_grid[2] is number of points) + _dx (float): horizontal size of the exit slit [m] + _dy (float): vertical size of the exit slit [m] + _re (float): photon energy range [eV] + _ne (int): number of photons vs. photon energy + _disp_rate (float): linear dispersion rate at the exit slit [eV/m] + _disp_pl (string): input character describing orientation of the dispersion plane ('h' or 'x' for horizontal, 'v' or 'y' for vertical, default) + _e0 (float): central relative photon energy for the resulting spectrum [eV] + _x0 (float): central horizontal position of the exit slit [m] + _y0 (float): central vertical position of the exit slit [m] + _nx_int (int): number of points vs. hor. position to use at integration within aperture + _ny_int (int): number of points vs. vert. position to use at integration within aperture + """ + + arSpec = array('d', [0]*_ne) + multFlux = 1.e+06 #Flux values multiplier, since input positions are in [m] whereas flux density in [ph/s/.1%bw/mm^2] + + #xc = 0.5*(_x_grid[0] + _x_grid[1]) + #yc = 0.5*(_y_grid[0] + _y_grid[1]) + dxHalf = 0.5*_dx + dyHalf = 0.5*_dy + xLim = [_x0-dxHalf, _x0+dxHalf, _nx_int] + yLim = [_y0-dyHalf, _y0+dyHalf, _ny_int] + + reHalf = 0.5*_re + posStart = (_e0 - reHalf)/_disp_rate + posEnd = (_e0 + reHalf)/_disp_rate + posStep = (posEnd - posStart)/(_ne - 1) + pos = posStart + + posIsY = False + if _disp_pl in ['v', 'V', 'y', 'Y']: posIsY = True + + for ie in range(_ne): + + if posIsY: + yLim[0] = pos - dyHalf + yLim[1] = pos + dyHalf + else: + xLim[0] = pos - dxHalf + xLim[1] = pos + dxHalf + + arSpec[ie] = uti_math.integ_ar_2d(_ar, 1, _x_grid, _y_grid, _x_lim=xLim, _y_lim=yLim)*multFlux + pos += posStep + + return arSpec + \ No newline at end of file diff --git a/env/python/srwpy/srwlib.py b/env/python/srwpy/srwlib.py index f1648487..bed9bb19 100644 --- a/env/python/srwpy/srwlib.py +++ b/env/python/srwpy/srwlib.py @@ -593,6 +593,31 @@ def add(self, _mag, _xc=None, _yc=None, _zc=None, _vx=None, _vy=None, _vz=None, self.arVy.append(_vy) self.arVz.append(_vz) self.arAng.append(_ang) + + def combine(self, _magCnt): #OC07092024 + if(_magCnt is None): return + if(not isinstance(_magCnt, SRWLMagFldC)): raise Exception("Magnetic field container is expected.") + + nOldElem = len(self.arMagFld) + nNewElem = len(_magCnt.arMagFld) + for i in range(nNewElem): + self.arMagFld.append(_magCnt.arMagFld[i]) + self.arXc.append(_magCnt.arXc[i]) + self.arYc.append(_magCnt.arYc[i]) + self.arZc.append(_magCnt.arZc[i]) + + if hasattr(_magCnt, 'arVx'): + if not hasattr(self, 'arVx'): self.arVx = array('d', [0]*nOldElem) + self.arVx.extend(_magCnt.arVx) + if hasattr(_magCnt, 'arVy'): + if not hasattr(self, 'arVy'): self.arVy = array('d', [0]*nOldElem) + self.arVy.extend(_magCnt.arVy) + if hasattr(_magCnt, 'arVz'): + if not hasattr(self, 'arVz'): self.arVz = array('d', [1]*nOldElem) + self.arVz.extend(_magCnt.arVz) + if hasattr(_magCnt, 'arAng'): + if not hasattr(self, 'arAng'): self.arAng = array('d', [0]*nOldElem) + self.arAng.extend(_magCnt.arAng) #**************************************************************************** class SRWLPrtTrj(object): diff --git a/env/python/srwpy/uti_math_eigen.py b/env/python/srwpy/uti_math_eigen.py index 14130bf5..2d486445 100644 --- a/env/python/srwpy/uti_math_eigen.py +++ b/env/python/srwpy/uti_math_eigen.py @@ -43,8 +43,12 @@ except: print("UtiMathEigen WARNING: sklearn unavailable; to install: pip install -U scikit-learn") import math -import numpy as np +#import numpy as np #OC23062024 (commented-out) +try: #OC23062024 + import numpy as np +except: + pass #class Linear_alg: class UtiMathEigen: #OC12052021