parameters: (name, x)
return_statement: return chi2
docstring: %SimTrans calculates transmission spectrum from two crystals
      %   lam - array containing wavelengths to calc over
      %   hkl - contains all Nref hkl's that calculation is performed for
      %   bgd - array containing coefficients of polynomial for background
      %   sf - overall scale factor
      %   pktype - 1 = gauss; 2 = lorentz; ...
      %   UB1 - UB matrix for first crystal
      %   setang1 - setting angles for first crystal (deviations from ideal UB
      %               location).
      %   pkpars1 - position(lambda), position(d-spacing), position(ttheta), width, intensity for each Nref reflections
      %   UB2,setang2,pkpars2 - as above, for second crystal
      %
      % M. Guthrie 21st Jan 2014
      %
      % calculate background profile
      % determine number of coeffs bgd
docstring_summary: %SimTrans calculates transmission spectrum from two crystals
      %   lam - array containing wavelengths to calc over
      %   hkl - contains all Nref hkl's that calculation is performed for
      %   bgd - array containing coefficients of polynomial for background
      %   sf - overall scale factor
      %   pktype - 1 = gauss; 2 = lorentz; ...
      %   UB1 - UB matrix for first crystal
      %   setang1 - setting angles for first crystal (deviations from ideal UB
      %               location).
      %   pkpars1 - position(lambda), position(d-spacing), position(ttheta), width, intensity for each Nref reflections
      %   UB2,setang2,pkpars2 - as above, for second crystal
      %
      % M. Guthrie 21st Jan 2014
      %
      % calculate background profile
      % determine number of coeffs bgd
docstring_tokens: %SimTrans calculates transmission spectrum from two crystals % lam - array containing wavelengths to calc over % hkl - contains all Nref hkl s that calculation is performed for % bgd - array containing coefficients of polynomial for background % sf - overall scale factor % pktype - 1 = gauss ; 2 = lorentz ; ... % UB1 - UB matrix for first crystal % setang1 - setting angles for first crystal ( deviations from ideal UB % location ) . % pkpars1 - position ( lambda ) position ( d - spacing ) position ( ttheta ) width intensity for each Nref reflections % UB2 setang2 pkpars2 - as above for second crystal % % M . Guthrie 21st Jan 2014 % % calculate background profile % determine number of coeffs bgd
function: def SimTransOutput3(name, x):
      '''
      %SimTrans calculates transmission spectrum from two crystals
      %   lam - array containing wavelengths to calc over
      %   hkl - contains all Nref hkl's that calculation is performed for
      %   bgd - array containing coefficients of polynomial for background
      %   sf - overall scale factor
      %   pktype - 1 = gauss; 2 = lorentz; ...
      %   UB1 - UB matrix for first crystal
      %   setang1 - setting angles for first crystal (deviations from ideal UB
      %               location).
      %   pkpars1 - position(lambda), position(d-spacing), position(ttheta), width, intensity for each Nref reflections
      %   UB2,setang2,pkpars2 - as above, for second crystal
      %
      % M. Guthrie 21st Jan 2014
      %
      % calculate background profile
      % determine number of coeffs bgd
      '''
      global hkl1, hkl2
      global UB1, pkcalcint1
      global UB2, pkcalcint2
      global pktype
      global lam, y, e, TOF
      global L1
      global ttot
      global fxsamediam
      global neqv1, eqvlab1, neqv2, eqvlab2
      global difa, function_verbose
      global figure_name_attenuation, run_number
  
      nref1 = hkl1.shape[0]  # % number of reflections to integrate over
      nref2 = hkl2.shape[0]  # % number of reflections to integrate over
      # % returns array with same dim as input labelling equivs
      eqvlab1, neqv1 = findeqvs(hkl1)
      eqvlab2, neqv2 = findeqvs(hkl2)
  
      setang1 = x[0:3]
      pkmult1 = x[3:4 + neqv1 - 1]
      setang2 = x[4 + neqv1 - 1:6 + neqv1]
      sf = x[neqv1 + neqv2 + 7 - 1]
      pkwid1 = x[neqv1 + neqv2 + 7]
      bgd = x[neqv1 + neqv2 + 8:neqv1 + neqv2 + 9 + 2 - 1]
      pkwid2 = x[neqv1 + neqv2 + 10]
      # % if diamond intensities the same, allow single scale f
      relsf = x[neqv1 + neqv2 + 11]
      if fxsamediam == 1:
          x[6 + neqv1:7 + neqv1 + neqv2 - 1] = x[3:4 + neqv2 - 1] * relsf
          pkmult2 = x[6 + neqv1:7 + neqv1 + neqv2 - 1]
      else:
          pkmult2 = x[6 + neqv1:7 + neqv1 + neqv2 - 1]
      delam = x[neqv1 + neqv2 + 12]
      L2 = x[neqv1 + neqv2 + 13]
  
      shftlam = 0.0039558 * TOF / (L1 + L2) + difa * (TOF ** 2)
      # number of lambda points to calculate over
      npt = shftlam.shape[0]
      # calculate information for peaks for crystal 1 using hkl,UB1, setang,
      # pkpos
      a, b, c = pkposcalc(hkl1, UB1, setang1)
      pkpars1 = np.column_stack((a, b, c))
      # calculate information for peaks for crystal 2 using hkl,UB1, setang,
      # pkpos
      a, b, c = pkposcalc(hkl2, UB2, setang2)
      pkpars2 = np.column_stack((a, b, c))
  
      # generate nptx,nco array containing, x^0,x^1,x^2,...x^nco for
      # all nonzero background coefficients
      bgdco = np.where(bgd != 0)[0]
      nco = bgdco.shape[0]
      nonzerobgd = np.zeros(nco)
  
      X = np.ones(shape=(nco, npt))
      for i in range(nco):
          X[i, :] = shftlam ** (bgd[bgdco[i]] - 1)
          nonzerobgd[i] = bgd[bgdco[i]]
  
      # calculate background profile by multiplying this with coefficients
      # themselves
      bgdprof = nonzerobgd.dot(X)
      # bgdprof = np.outer(nonzerobgd, X)
      # print bgdprof
      # bgdprof = bgdprof[0, :]
      # calculate peaks for crystal 1
  
      t1 = np.zeros(npt)  # initialise array containing profile
      for i in range(nref1):
          if pktype == 1:
              pkpars1[i][0] = pkpars1[i][0] * delam  # linear lambda shift
              sig = pkwid1 * pkpars1[i][0] + pkwid2 * (pkpars1[i][0] ** 2.)  # const del(lambda)/lambda
              extScl = pkpars1[i][0] ** 0  # lambda dependent extinction effect
              t1 = t1 - extScl * pkmult1[int(eqvlab1[i])] * pkcalcint1[i] * (
                  np.exp(-((shftlam - pkpars1[i][0]) ** 2.) / (2 * (sig ** 2))))
  
      # calculate peaks for crystal 2
      t2 = np.zeros(npt)  # initialise array containing profile
      for i in range(nref2):
          if pktype == 1:
              pkpars2[i][0] = pkpars2[i][0] * delam  # linear lambda shift
              sig = pkwid1 * pkpars2[i][0] + pkwid2 * (pkpars2[i][0] ** 2.)  # const del(lambda)/lambda
              extScl = pkpars2[i][0] ** 0  # lambda dependent extinction effect
              t2 = t2 - extScl * pkmult2[int(eqvlab2[i])] * pkcalcint2[i] * (
                  np.exp(-(shftlam - pkpars2[i][0]) ** 2. / (2 * (sig ** 2))))
  
      # calculate final profile
      ttot = (bgdprof + sf * t1) * (bgdprof + sf * t2)
      # t1 = 1.0;
      # t2 = 1.0;
      # introduce weighting function and calc chi2...
      w = np.ones(len(shftlam))  # equal weighting everywhere
      # i1 = np.where(shftlam > 2.15)[0][0]
      # j1 = np.where(shftlam > 2.65)[0][0]
      # w[i1:j1] = 5 #extra weighting in region of first peaks
      # i1 = find(lam>1.68,1,'first');
      # j1 = find(lam>2.05,1,'first');
      # w(i1:j1)=5; %extra weighting but not too much
  
      resid = (y - ttot) * w
      chi2 = np.sum(resid ** 2. / (2 * e ** 2)) / npt
  
      output = 1
      if output == 1:
          diam1trans = sf * t1 + bgdprof
          diam2trans = sf * t2 + bgdprof
          out = np.column_stack((shftlam, diam1trans, diam2trans, ttot, y))
          np.savetxt(name, out, delimiter=',')
  
      figure_name_attenuation = 'Attenuation ' + run_number
  
      plt.figure(figure_name_attenuation)
      plt.plot(shftlam, ttot, 'r', label='Total att.')
      plt.plot(shftlam, diam1trans, 'k', label='Diam 1 att.')
      plt.plot(shftlam, diam2trans, 'b', label='Diam 2 att.')
      plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
                 ncol=3, mode="expand", borderaxespad=0.)
      plt.xlabel('Wavelength (A)')
      plt.ylabel('Transmission')
      plt.grid()
      for i in range(len(pkpars1)):
          plt.arrow(pkpars1[i, 0] * delam, 1.1, 0.0, 0.025,
                    fc="k", ec="k", head_width=0, head_length=0)
      for i in range(len(pkpars2)):
          plt.arrow(pkpars2[i, 0] * delam, 1.15, 0.0, 0.025,
                    fc="k", ec="k", head_width=0, head_length=0)
      plt.xlim(1, 2.7)
      plt.ylim(0, 1.2)
      plt.show()
  
      return chi2
function_tokens: def SimTransOutput3 ( name , x ) : global hkl1 , hkl2 global UB1 , pkcalcint1 global UB2 , pkcalcint2 global pktype global lam , y , e , TOF global L1 global ttot global fxsamediam global neqv1 , eqvlab1 , neqv2 , eqvlab2 global difa , function_verbose global figure_name_attenuation , run_number nref1 = hkl1 . shape [ 0 ] # % number of reflections to integrate over nref2 = hkl2 . shape [ 0 ] # % number of reflections to integrate over # % returns array with same dim as input labelling equivs eqvlab1 , neqv1 = findeqvs ( hkl1 ) eqvlab2 , neqv2 = findeqvs ( hkl2 ) setang1 = x [ 0 : 3 ] pkmult1 = x [ 3 : 4 + neqv1 - 1 ] setang2 = x [ 4 + neqv1 - 1 : 6 + neqv1 ] sf = x [ neqv1 + neqv2 + 7 - 1 ] pkwid1 = x [ neqv1 + neqv2 + 7 ] bgd = x [ neqv1 + neqv2 + 8 : neqv1 + neqv2 + 9 + 2 - 1 ] pkwid2 = x [ neqv1 + neqv2 + 10 ] # % if diamond intensities the same, allow single scale f relsf = x [ neqv1 + neqv2 + 11 ] if fxsamediam == 1 : x [ 6 + neqv1 : 7 + neqv1 + neqv2 - 1 ] = x [ 3 : 4 + neqv2 - 1 ] * relsf pkmult2 = x [ 6 + neqv1 : 7 + neqv1 + neqv2 - 1 ] else : pkmult2 = x [ 6 + neqv1 : 7 + neqv1 + neqv2 - 1 ] delam = x [ neqv1 + neqv2 + 12 ] L2 = x [ neqv1 + neqv2 + 13 ] shftlam = 0.0039558 * TOF / ( L1 + L2 ) + difa * ( TOF ** 2 ) # number of lambda points to calculate over npt = shftlam . shape [ 0 ] # calculate information for peaks for crystal 1 using hkl,UB1, setang, # pkpos a , b , c = pkposcalc ( hkl1 , UB1 , setang1 ) pkpars1 = np . column_stack ( ( a , b , c ) ) # calculate information for peaks for crystal 2 using hkl,UB1, setang, # pkpos a , b , c = pkposcalc ( hkl2 , UB2 , setang2 ) pkpars2 = np . column_stack ( ( a , b , c ) ) # generate nptx,nco array containing, x^0,x^1,x^2,...x^nco for # all nonzero background coefficients bgdco = np . where ( bgd != 0 ) [ 0 ] nco = bgdco . shape [ 0 ] nonzerobgd = np . zeros ( nco ) X = np . ones ( shape = ( nco , npt ) ) for i in range ( nco ) : X [ i , : ] = shftlam ** ( bgd [ bgdco [ i ] ] - 1 ) nonzerobgd [ i ] = bgd [ bgdco [ i ] ] # calculate background profile by multiplying this with coefficients # themselves bgdprof = nonzerobgd . dot ( X ) # bgdprof = np.outer(nonzerobgd, X) # print bgdprof # bgdprof = bgdprof[0, :] # calculate peaks for crystal 1 t1 = np . zeros ( npt ) # initialise array containing profile for i in range ( nref1 ) : if pktype == 1 : pkpars1 [ i ] [ 0 ] = pkpars1 [ i ] [ 0 ] * delam # linear lambda shift sig = pkwid1 * pkpars1 [ i ] [ 0 ] + pkwid2 * ( pkpars1 [ i ] [ 0 ] ** 2. ) # const del(lambda)/lambda extScl = pkpars1 [ i ] [ 0 ] ** 0 # lambda dependent extinction effect t1 = t1 - extScl * pkmult1 [ int ( eqvlab1 [ i ] ) ] * pkcalcint1 [ i ] * ( np . exp ( - ( ( shftlam - pkpars1 [ i ] [ 0 ] ) ** 2. ) / ( 2 * ( sig ** 2 ) ) ) ) # calculate peaks for crystal 2 t2 = np . zeros ( npt ) # initialise array containing profile for i in range ( nref2 ) : if pktype == 1 : pkpars2 [ i ] [ 0 ] = pkpars2 [ i ] [ 0 ] * delam # linear lambda shift sig = pkwid1 * pkpars2 [ i ] [ 0 ] + pkwid2 * ( pkpars2 [ i ] [ 0 ] ** 2. ) # const del(lambda)/lambda extScl = pkpars2 [ i ] [ 0 ] ** 0 # lambda dependent extinction effect t2 = t2 - extScl * pkmult2 [ int ( eqvlab2 [ i ] ) ] * pkcalcint2 [ i ] * ( np . exp ( - ( shftlam - pkpars2 [ i ] [ 0 ] ) ** 2. / ( 2 * ( sig ** 2 ) ) ) ) # calculate final profile ttot = ( bgdprof + sf * t1 ) * ( bgdprof + sf * t2 ) # t1 = 1.0; # t2 = 1.0; # introduce weighting function and calc chi2... w = np . ones ( len ( shftlam ) ) # equal weighting everywhere # i1 = np.where(shftlam > 2.15)[0][0] # j1 = np.where(shftlam > 2.65)[0][0] # w[i1:j1] = 5 #extra weighting in region of first peaks # i1 = find(lam>1.68,1,'first'); # j1 = find(lam>2.05,1,'first'); # w(i1:j1)=5; %extra weighting but not too much resid = ( y - ttot ) * w chi2 = np . sum ( resid ** 2. / ( 2 * e ** 2 ) ) / npt output = 1 if output == 1 : diam1trans = sf * t1 + bgdprof diam2trans = sf * t2 + bgdprof out = np . column_stack ( ( shftlam , diam1trans , diam2trans , ttot , y ) ) np . savetxt ( name , out , delimiter = ',' ) figure_name_attenuation = 'Attenuation ' + run_number plt . figure ( figure_name_attenuation ) plt . plot ( shftlam , ttot , 'r' , label = 'Total att.' ) plt . plot ( shftlam , diam1trans , 'k' , label = 'Diam 1 att.' ) plt . plot ( shftlam , diam2trans , 'b' , label = 'Diam 2 att.' ) plt . legend ( bbox_to_anchor = ( 0. , 1.02 , 1. , .102 ) , loc = 3 , ncol = 3 , mode = "expand" , borderaxespad = 0. ) plt . xlabel ( 'Wavelength (A)' ) plt . ylabel ( 'Transmission' ) plt . grid ( ) for i in range ( len ( pkpars1 ) ) : plt . arrow ( pkpars1 [ i , 0 ] * delam , 1.1 , 0.0 , 0.025 , fc = "k" , ec = "k" , head_width = 0 , head_length = 0 ) for i in range ( len ( pkpars2 ) ) : plt . arrow ( pkpars2 [ i , 0 ] * delam , 1.15 , 0.0 , 0.025 , fc = "k" , ec = "k" , head_width = 0 , head_length = 0 ) plt . xlim ( 1 , 2.7 ) plt . ylim ( 0 , 1.2 ) plt . show ( ) return chi2