[G8b_run] Polarization modification code

Michael Dugger dugger at jlab.org
Wed Nov 9 18:29:09 EST 2011


Hi,

I have some code that I call "polValNew code" that takes
the polarization from Ken's tables and modifies it to a
value named polValNew.

Dependencies:

The polValNew code calls a function named polFix.

The polFix function calls a function named polFix0.

The polFix0 function requires a common block vector
named ldPar(6). Note: The vector is of length 6, but
only 4 values are used in the code.

The common block vector ldPar has elements that depend
on the run conditions: Manual/Auto, coherent edge,
polarization plane.

Even for those that do not use FORTRAN, the code is simple
and should be easy to read.

All of the required information can be found below:

c--> polValNew code:
###################
       eDiff = (edgeVec(i) - egVec(j))/1000

c--> Note: edgeVec(i) is the coherent edge in MeV
c--> Note: egVec(j) is the incident photon energy in MeV

       polValNew = polVal

c--> Note: polVal is the polarization value from Ken's table
c--> Note: I get the best results for diffFix = 0.135
c--> Note: Use plane = 1 for PARA and plane = 2 for PERP

       if (polFix(eDiff,plane,diffFix).gt.0) then
          polValNew = polVal/polFix(eDiff,plane,diffFix)
       endif

###################

c--> polFix.f code:
###################
       REAL function polFix(xE,type,diffFix)
       REAL xE,norm,diffFix
       INTEGER type
       polFix = polFix0(xE,type)
       norm = polFix0(diffFix,type)
       polFix = polFix/norm
       END
###################

c--> polFix0.f code:
###################
       REAL function polFix0(xE,type)
       REAL ldPar(6),fun,xE,yVal,m,b,par(3)
       REAL par1,par2,par3,tmpFun
       INTEGER i,j,type,region
       COMMON /ldBlock/ ldPar

       region = -1
       if (xE.le.0.025) region = 0
       if (xE.gt.0.025.and.xE.le.0.050) region = 1
       if (xE.gt.0.050.and.xE.le.0.075) region = 2
       if (xE.gt.0.075.and.xE.le.0.100) region = 3
       if (xE.gt.0.100.and.xE.le.0.125) region = 4
       if (xE.gt.0.125.and.xE.le.0.150) region = 5
       if (xE.gt.0.150.and.xE.le.0.175) region = 6
       if (xE.gt.0.175.and.xE.le.0.200) region = 7

       m = 1.0/0.025
       b = -0.5
       yVal = m*xE + b

       if (type.eq.1) then
          par(1) = ldPar(1)
          par(2) = ldPar(2)
       else
          par(1) = ldPar(4)
          par(2) = ldPar(5)
       endif

       par1 = par(1)
       par2 = par(2)

       tmpFun = 1.0
       if (region.ge.1) then
          do i=1,region
             tmpFun = tmpFun*(par1 + par2*yVal)
             yVal = yVal - 1.0
          enddo
       endif
       polFix0 = tmpFun

       if (tmpFun.eq.0) polFix0 = 1.0

       END
###################

ldPar values:

For 1.3 Manual:
ldPar(1) = 1.004337
ldPar(2) = -.4986358E-02
ldPar(4) = 1.009058
ldPar(5) = -.5687501E-02

For 1.5 Manual:
ldPar(1) = 1.012704
ldPar(2) = -.3850369E-02
ldPar(4) = 1.002215
ldPar(5) = -.2256782E-02

For 1.7 Manual:
ldPar(1) = 0.9888860
ldPar(2) = -.2256782E-02
ldPar(4) = 0.9924032
ldPar(5) = -.3019766E-02

For 2.1 Manual:
ldPar(1) = 1.004469
ldPar(2) = -.2188895E-01
ldPar(4) = 1.067167
ldPar(5) = -.1876592E-01

For 1.7 Auto:
ldPar(1) = 0.9844326
ldPar(2) = 0.7438295E-03
ldPar(4) = 1.032343
ldPar(5) = -.1255454E-01

For 1.9 Auto:
ldPar(1) = 1.008734
ldPar(2) = -.1557868E-02
ldPar(4) = 0.9865551
ldPar(5) = 0.2070320E-02

###################

Please let me know if there are any questions or if it appears that I have 
left something out.

Take care,
Michael



More information about the G8b_run mailing list