Context Navigation

source: trunk/GSASIIstrMath.py @ 1888

Last change on this file since 1888 was 1888, checked in by vondreele, 10 years ago

Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 146.1 KB

Rev	Line
[1299]	1	# -- coding: utf-8 --
	2	'''
	3	GSASIIstrMath - structure math routines
	4	-----------------------------------------
	5	'''
	6	########### SVN repository information ###################
	7	# $Date: 2015-06-15 13:56:19 +0000 (Mon, 15 Jun 2015) $
	8	# $Author: vondreele $
	9	# $Revision: 1888 $
	10	# $URL: trunk/GSASIIstrMath.py $
	11	# $Id: GSASIIstrMath.py 1888 2015-06-15 13:56:19Z vondreele $
	12	########### SVN repository information ###################
	13	import time
	14	import copy
	15	import numpy as np
	16	import numpy.ma as ma
	17	import numpy.linalg as nl
	18	import scipy.optimize as so
	19	import scipy.stats as st
	20	import GSASIIpath
	21	GSASIIpath.SetVersionNumber("$Revision: 1888 $")
	22	import GSASIIElem as G2el
	23	import GSASIIlattice as G2lat
	24	import GSASIIspc as G2spc
	25	import GSASIIpwd as G2pwd
	26	import GSASIImapvars as G2mv
	27	import GSASIImath as G2mth
[1813]	28	import GSASIIobj as G2obj
[1299]	29
	30	sind = lambda x: np.sin(x*np.pi/180.)
	31	cosd = lambda x: np.cos(x*np.pi/180.)
	32	tand = lambda x: np.tan(x*np.pi/180.)
	33	asind = lambda x: 180.*np.arcsin(x)/np.pi
	34	acosd = lambda x: 180.*np.arccos(x)/np.pi
	35	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
	36
[1884]	37	ateln2 = 8.0*np.log(2.0)
[1630]	38	twopi = 2.0*np.pi
	39	twopisq = 2.0np.pi*2
[1299]	40
	41	################################################################################
	42	##### Rigid Body Models
	43	################################################################################
	44
	45	def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
	46	''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
	47	current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
	48	'''
	49	atxIds = ['Ax:','Ay:','Az:']
	50	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
	51	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
	52	if not RBIds['Vector'] and not RBIds['Residue']:
	53	return
	54	VRBIds = RBIds['Vector']
	55	RRBIds = RBIds['Residue']
	56	if Update:
	57	RBData = rigidbodyDict
	58	else:
	59	RBData = copy.deepcopy(rigidbodyDict) # don't mess with original!
	60	if RBIds['Vector']: # first update the vector magnitudes
	61	VRBData = RBData['Vector']
	62	for i,rbId in enumerate(VRBIds):
	63	if VRBData[rbId]['useCount']:
	64	for j in range(len(VRBData[rbId]['VectMag'])):
	65	name = '::RBV;'+str(j)+':'+str(i)
	66	VRBData[rbId]['VectMag'][j] = parmDict[name]
	67	for phase in Phases:
	68	Phase = Phases[phase]
	69	General = Phase['General']
[1604]	70	cx,ct,cs,cia = General['AtomPtrs']
[1299]	71	cell = General['Cell'][1:7]
	72	Amat,Bmat = G2lat.cell2AB(cell)
[1604]	73	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
[1299]	74	pfx = str(Phase['pId'])+'::'
	75	if Update:
	76	RBModels = Phase['RBModels']
	77	else:
	78	RBModels = copy.deepcopy(Phase['RBModels']) # again don't mess with original!
	79	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
	80	jrb = VRBIds.index(RBObj['RBId'])
	81	rbsx = str(irb)+':'+str(jrb)
	82	for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
	83	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
	84	for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
	85	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
	86	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
	87	TLS = RBObj['ThermalMotion']
	88	if 'T' in TLS[0]:
	89	for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
	90	TLS[1][i] = parmDict[pfx+pt+rbsx]
	91	if 'L' in TLS[0]:
	92	for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
	93	TLS[1][i+6] = parmDict[pfx+pt+rbsx]
	94	if 'S' in TLS[0]:
	95	for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
	96	TLS[1][i+12] = parmDict[pfx+pt+rbsx]
	97	if 'U' in TLS[0]:
	98	TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
	99	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
	100	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
	101	for i,x in enumerate(XYZ):
	102	atId = RBObj['Ids'][i]
	103	for j in [0,1,2]:
	104	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
	105	if UIJ[i][0] == 'A':
	106	for j in range(6):
	107	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
	108	elif UIJ[i][0] == 'I':
	109	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
	110
	111	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
	112	jrb = RRBIds.index(RBObj['RBId'])
	113	rbsx = str(irb)+':'+str(jrb)
	114	for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
	115	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
	116	for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
	117	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
	118	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
	119	TLS = RBObj['ThermalMotion']
	120	if 'T' in TLS[0]:
	121	for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
	122	RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
	123	if 'L' in TLS[0]:
	124	for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
	125	RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
	126	if 'S' in TLS[0]:
	127	for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
	128	RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
	129	if 'U' in TLS[0]:
	130	RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
	131	for itors,tors in enumerate(RBObj['Torsions']):
	132	tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
	133	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
	134	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
	135	for i,x in enumerate(XYZ):
	136	atId = RBObj['Ids'][i]
	137	for j in [0,1,2]:
	138	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
	139	if UIJ[i][0] == 'A':
	140	for j in range(6):
	141	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
	142	elif UIJ[i][0] == 'I':
	143	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
	144
	145	def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
	146	'Needs a doc string'
	147	atxIds = ['dAx:','dAy:','dAz:']
	148	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
	149	TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
	150	LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
	151	SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
	152	PIds = ['Px:','Py:','Pz:']
	153	OIds = ['Oa:','Oi:','Oj:','Ok:']
	154	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
	155	if not RBIds['Vector'] and not RBIds['Residue']:
	156	return
	157	VRBIds = RBIds['Vector']
	158	RRBIds = RBIds['Residue']
	159	RBData = rigidbodyDict
	160	for item in parmDict:
	161	if 'RB' in item:
	162	dFdvDict[item] = 0. #NB: this is a vector which is no. refl. long & must be filled!
	163	General = Phase['General']
[1604]	164	cx,ct,cs,cia = General['AtomPtrs']
[1299]	165	cell = General['Cell'][1:7]
	166	Amat,Bmat = G2lat.cell2AB(cell)
	167	rpd = np.pi/180.
	168	rpd2 = rpd**2
	169	g = nl.inv(np.inner(Bmat,Bmat))
	170	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
	171	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
[1604]	172	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
[1299]	173	pfx = str(Phase['pId'])+'::'
	174	RBModels = Phase['RBModels']
	175	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
	176	VModel = RBData['Vector'][RBObj['RBId']]
	177	Q = RBObj['Orient'][0]
	178	Pos = RBObj['Orig'][0]
	179	jrb = VRBIds.index(RBObj['RBId'])
	180	rbsx = str(irb)+':'+str(jrb)
	181	dXdv = []
	182	for iv in range(len(VModel['VectMag'])):
	183	dCdv = []
	184	for vec in VModel['rbVect'][iv]:
	185	dCdv.append(G2mth.prodQVQ(Q,vec))
	186	dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
	187	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
	188	for ia,atId in enumerate(RBObj['Ids']):
	189	atNum = AtLookup[atId]
	190	dx = 0.00001
	191	for iv in range(len(VModel['VectMag'])):
	192	for ix in [0,1,2]:
	193	dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	194	for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
	195	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
	196	for iv in range(4):
	197	Q[iv] -= dx
	198	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
	199	Q[iv] += 2.*dx
	200	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
	201	Q[iv] -= dx
	202	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
	203	for ix in [0,1,2]:
	204	dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	205	X = G2mth.prodQVQ(Q,Cart[ia])
	206	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
	207	dFdu = G2lat.U6toUij(dFdu.T)
	208	dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
	209	dFdu = G2lat.UijtoU6(dFdu)
	210	atNum = AtLookup[atId]
	211	if 'T' in RBObj['ThermalMotion'][0]:
	212	for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
	213	dFdvDict[pfx+name+rbsx] += dFdu[i]
	214	if 'L' in RBObj['ThermalMotion'][0]:
	215	dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
	216	dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
	217	dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
	218	dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
	219	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
	220	dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2(-dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
	221	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
	222	dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2(-dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
	223	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
	224	if 'S' in RBObj['ThermalMotion'][0]:
	225	dFdvDict[pfx+'RBVS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
	226	dFdvDict[pfx+'RBVS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
	227	dFdvDict[pfx+'RBVS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
	228	dFdvDict[pfx+'RBVS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
	229	dFdvDict[pfx+'RBVS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
	230	dFdvDict[pfx+'RBVS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
	231	dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
	232	dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
	233	if 'U' in RBObj['ThermalMotion'][0]:
	234	dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
	235
	236
	237	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
	238	Q = RBObj['Orient'][0]
	239	Pos = RBObj['Orig'][0]
	240	jrb = RRBIds.index(RBObj['RBId'])
	241	torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
	242	rbsx = str(irb)+':'+str(jrb)
	243	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
	244	for itors,tors in enumerate(RBObj['Torsions']): #derivative error?
	245	tname = pfx+'RBRTr;'+str(itors)+':'+rbsx
	246	orId,pvId = torData[itors][:2]
	247	pivotVec = Cart[orId]-Cart[pvId]
	248	QA = G2mth.AVdeg2Q(-0.001,pivotVec)
	249	QB = G2mth.AVdeg2Q(0.001,pivotVec)
	250	for ir in torData[itors][3]:
	251	atNum = AtLookup[RBObj['Ids'][ir]]
	252	rVec = Cart[ir]-Cart[pvId]
	253	dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
	254	dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
	255	for ix in [0,1,2]:
	256	dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	257	for ia,atId in enumerate(RBObj['Ids']):
	258	atNum = AtLookup[atId]
	259	dx = 0.00001
	260	for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
	261	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
	262	for iv in range(4):
	263	Q[iv] -= dx
	264	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
	265	Q[iv] += 2.*dx
	266	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
	267	Q[iv] -= dx
	268	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
	269	for ix in [0,1,2]:
	270	dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	271	X = G2mth.prodQVQ(Q,Cart[ia])
	272	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
	273	dFdu = G2lat.U6toUij(dFdu.T)
	274	dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
	275	dFdu = G2lat.UijtoU6(dFdu)
	276	atNum = AtLookup[atId]
	277	if 'T' in RBObj['ThermalMotion'][0]:
	278	for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
	279	dFdvDict[pfx+name+rbsx] += dFdu[i]
	280	if 'L' in RBObj['ThermalMotion'][0]:
	281	dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
	282	dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
	283	dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
	284	dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
	285	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
	286	dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2(dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
	287	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
	288	dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2(dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
	289	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
	290	if 'S' in RBObj['ThermalMotion'][0]:
	291	dFdvDict[pfx+'RBRS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
	292	dFdvDict[pfx+'RBRS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
	293	dFdvDict[pfx+'RBRS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
	294	dFdvDict[pfx+'RBRS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
	295	dFdvDict[pfx+'RBRS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
	296	dFdvDict[pfx+'RBRS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
	297	dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
	298	dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
	299	if 'U' in RBObj['ThermalMotion'][0]:
	300	dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
	301
	302	################################################################################
	303	##### Penalty & restraint functions
	304	################################################################################
	305
[1775]	306	def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
[1299]	307	'Needs a doc string'
	308	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	309	pNames = []
	310	pVals = []
	311	pWt = []
	312	negWt = {}
	313	pWsum = {}
	314	for phase in Phases:
	315	pId = Phases[phase]['pId']
	316	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
	317	General = Phases[phase]['General']
[1604]	318	cx,ct,cs,cia = General['AtomPtrs']
[1299]	319	textureData = General['SH Texture']
	320	SGData = General['SGData']
[1604]	321	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
[1299]	322	cell = General['Cell'][1:7]
	323	Amat,Bmat = G2lat.cell2AB(cell)
	324	if phase not in restraintDict:
	325	continue
	326	phaseRest = restraintDict[phase]
	327	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
	328	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
[1775]	329	['ChemComp','Sites'],['Texture','HKLs'],]
[1299]	330	for name,rest in names:
	331	pWsum[name] = 0.
	332	itemRest = phaseRest[name]
	333	if itemRest[rest] and itemRest['Use']:
	334	wt = itemRest['wtFactor']
	335	if name in ['Bond','Angle','Plane','Chiral']:
	336	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
	337	pNames.append(str(pId)+':'+name+':'+str(i))
	338	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	339	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
	340	if name == 'Bond':
	341	calc = G2mth.getRestDist(XYZ,Amat)
	342	elif name == 'Angle':
	343	calc = G2mth.getRestAngle(XYZ,Amat)
	344	elif name == 'Plane':
	345	calc = G2mth.getRestPlane(XYZ,Amat)
	346	elif name == 'Chiral':
	347	calc = G2mth.getRestChiral(XYZ,Amat)
	348	pVals.append(obs-calc)
	349	pWt.append(wt/esd**2)
	350	pWsum[name] += wt((obs-calc)/esd)*2
	351	elif name in ['Torsion','Rama']:
	352	coeffDict = itemRest['Coeff']
	353	for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
	354	pNames.append(str(pId)+':'+name+':'+str(i))
	355	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	356	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
	357	if name == 'Torsion':
	358	tor = G2mth.getRestTorsion(XYZ,Amat)
	359	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
	360	else:
	361	phi,psi = G2mth.getRestRama(XYZ,Amat)
	362	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
	363	pVals.append(restr)
	364	pWt.append(wt/esd**2)
	365	pWsum[name] += wt(restr/esd)*2
	366	elif name == 'ChemComp':
	367	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
	368	pNames.append(str(pId)+':'+name+':'+str(i))
	369	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
	370	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
	371	calc = np.sum(mulfracfactors)
	372	pVals.append(obs-calc)
	373	pWt.append(wt/esd**2)
	374	pWsum[name] += wt((obs-calc)/esd)*2
	375	elif name == 'Texture':
	376	SHkeys = textureData['SH Coeff'][1].keys()
	377	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
	378	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
	379	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
	380	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
	381	PH = np.array(hkl)
	382	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
	383	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
	384	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
	385	Z1 = -ma.masked_greater(Z,0.0)
	386	IndZ1 = np.array(ma.nonzero(Z1))
	387	for ind in IndZ1.T:
	388	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
	389	pVals.append(Z1[ind[0]][ind[1]])
	390	pWt.append(wt/esd1**2)
[1776]	391	pWsum[name] += wt(-Z1[ind[0]][ind[1]]/esd1)*2
[1299]	392	if ifesd2:
	393	Z2 = 1.-Z
	394	for ind in np.ndindex(grid,grid):
	395	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
	396	pVals.append(Z1[ind[0]][ind[1]])
	397	pWt.append(wt/esd2**2)
[1776]	398	pWsum[name] += wt(Z2/esd2)*2
[1775]	399
[1810]	400	for phase in Phases:
[1775]	401	name = 'SH-Pref.Ori.'
[1810]	402	pId = Phases[phase]['pId']
	403	General = Phases[phase]['General']
	404	SGData = General['SGData']
	405	cell = General['Cell'][1:7]
[1775]	406	pWsum[name] = 0.0
	407	for hist in Phases[phase]['Histograms']:
[1816]	408	if hist in Histograms and 'PWDR' in hist:
[1775]	409	hId = Histograms[hist]['hId']
	410	phfx = '%d:%d:'%(pId,hId)
	411	if calcControls[phfx+'poType'] == 'SH':
	412	toler = calcControls[phfx+'SHtoler']
	413	wt = 1./toler**2
[1793]	414	HKLs = np.array(calcControls[phfx+'SHhkl'])
[1777]	415	SHnames = calcControls[phfx+'SHnames']
	416	SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
[1775]	417	for i,PH in enumerate(HKLs):
	418	phi,beta = G2lat.CrsAng(PH,cell,SGData)
	419	SH3Coef = {}
	420	for item in SHcof:
	421	L,N = eval(item.strip('C'))
	422	SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]
	423	ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
	424	X = np.linspace(0,90.0,26)
	425	Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
	426	IndY = ma.nonzero(Y)
	427	for ind in IndY[0]:
	428	pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
	429	pVals.append(Y[ind])
	430	pWt.append(wt)
[1777]	431	pWsum[name] += wt(Y[ind])*2
[1299]	432	pWsum['PWLref'] = 0.
	433	for item in varyList:
	434	if 'PWLref' in item and parmDict[item] < 0.:
	435	pId = int(item.split(':')[0])
	436	if negWt[pId]:
	437	pNames.append(item)
	438	pVals.append(-parmDict[item])
	439	pWt.append(negWt[pId])
	440	pWsum['PWLref'] += negWt[pId](-parmDict[item])*2
	441	pVals = np.array(pVals)
	442	pWt = np.array(pWt) #should this be np.sqrt?
	443	return pNames,pVals,pWt,pWsum
	444
[1775]	445	def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
[1299]	446	'Needs a doc string'
	447	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	448	pDerv = np.zeros((len(varyList),len(pVal)))
	449	for phase in Phases:
	450	# if phase not in restraintDict:
	451	# continue
	452	pId = Phases[phase]['pId']
	453	General = Phases[phase]['General']
[1604]	454	cx,ct,cs,cia = General['AtomPtrs']
[1299]	455	SGData = General['SGData']
[1604]	456	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
[1299]	457	cell = General['Cell'][1:7]
	458	Amat,Bmat = G2lat.cell2AB(cell)
	459	textureData = General['SH Texture']
	460
	461	SHkeys = textureData['SH Coeff'][1].keys()
	462	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
	463	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
	464	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
	465	sam = SamSym[textureData['Model']]
	466	phaseRest = restraintDict.get(phase,{})
	467	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
	468	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
	469	'ChemComp':'Sites','Texture':'HKLs'}
	470	lasthkl = np.array([0,0,0])
	471	for ip,pName in enumerate(pNames):
	472	pnames = pName.split(':')
	473	if pId == int(pnames[0]):
	474	name = pnames[1]
	475	if 'PWL' in pName:
	476	pDerv[varyList.index(pName)][ip] += 1.
	477	continue
[1775]	478	elif 'SH-' in pName:
	479	continue
[1299]	480	id = int(pnames[2])
	481	itemRest = phaseRest[name]
	482	if name in ['Bond','Angle','Plane','Chiral']:
	483	indx,ops,obs,esd = itemRest[names[name]][id]
	484	dNames = []
	485	for ind in indx:
	486	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
	487	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	488	if name == 'Bond':
	489	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
	490	elif name == 'Angle':
	491	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
	492	elif name == 'Plane':
	493	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
	494	elif name == 'Chiral':
	495	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
	496	elif name in ['Torsion','Rama']:
	497	coffDict = itemRest['Coeff']
	498	indx,ops,cofName,esd = itemRest[names[name]][id]
	499	dNames = []
	500	for ind in indx:
	501	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
	502	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	503	if name == 'Torsion':
	504	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
	505	else:
	506	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
	507	elif name == 'ChemComp':
	508	indx,factors,obs,esd = itemRest[names[name]][id]
	509	dNames = []
	510	for ind in indx:
	511	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
	512	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
	513	deriv = mul*factors
	514	elif 'Texture' in name:
	515	deriv = []
	516	dNames = []
	517	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
	518	hkl = np.array(hkl)
	519	if np.any(lasthkl-hkl):
	520	PH = np.array(hkl)
	521	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
	522	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
	523	lasthkl = copy.copy(hkl)
	524	if 'unit' in name:
	525	pass
	526	else:
	527	gam = float(pnames[3])
	528	psi = float(pnames[4])
	529	for SHname in ODFln:
	530	l,m,n = eval(SHname[1:])
	531	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
	532	dNames += [str(pId)+'::'+SHname]
	533	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
	534	for dName,drv in zip(dNames,deriv):
	535	try:
	536	ind = varyList.index(dName)
	537	pDerv[ind][ip] += drv
	538	except ValueError:
	539	pass
[1775]	540
[1777]	541	lasthkl = np.array([0,0,0])
	542	for ip,pName in enumerate(pNames):
	543	deriv = []
	544	dNames = []
	545	pnames = pName.split(':')
	546	if 'SH-' in pName and pId == int(pnames[0]):
	547	hId = int(pnames[1])
	548	phfx = '%d:%d:'%(pId,hId)
	549	psi = float(pnames[4])
	550	HKLs = calcControls[phfx+'SHhkl']
	551	SHnames = calcControls[phfx+'SHnames']
	552	SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
	553	hkl = np.array(HKLs[int(pnames[3])])
	554	if np.any(lasthkl-hkl):
	555	PH = np.array(hkl)
	556	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
	557	SH3Coef = {}
	558	for item in SHcof:
	559	L,N = eval(item.strip('C'))
	560	SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]
	561	ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
	562	lasthkl = copy.copy(hkl)
	563	for SHname in SHnames:
	564	l,n = eval(SHname[1:])
	565	SH3name = 'C%d,0,%d'%(l,n)
	566	Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
	567	dNames += [phfx+SHname]
	568	deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
	569	for dName,drv in zip(dNames,deriv):
	570	try:
	571	ind = varyList.index(dName)
	572	pDerv[ind][ip] += drv
	573	except ValueError:
	574	pass
[1299]	575	return pDerv
	576
	577	################################################################################
	578	##### Function & derivative calculations
	579	################################################################################
	580
	581	def GetAtomFXU(pfx,calcControls,parmDict):
	582	'Needs a doc string'
	583	Natoms = calcControls['Natoms'][pfx]
	584	Tdata = Natoms*[' ',]
	585	Mdata = np.zeros(Natoms)
	586	IAdata = Natoms*[' ',]
	587	Fdata = np.zeros(Natoms)
	588	FFdata = []
	589	BLdata = []
	590	Xdata = np.zeros((3,Natoms))
	591	dXdata = np.zeros((3,Natoms))
	592	Uisodata = np.zeros(Natoms)
	593	Uijdata = np.zeros((6,Natoms))
	594	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
	595	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
	596	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
	597	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
	598	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
	599	for iatm in range(Natoms):
	600	for key in keys:
	601	parm = pfx+key+str(iatm)
	602	if parm in parmDict:
	603	keys[key][iatm] = parmDict[parm]
[1308]	604	Fdata = np.where(Fdata,Fdata,1.e-8) #avoid divide by zero in derivative calc.?
[1299]	605	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
	606
[1618]	607	def GetAtomSSFXU(pfx,calcControls,parmDict):
	608	'Needs a doc string'
	609	Natoms = calcControls['Natoms'][pfx]
[1625]	610	maxSSwave = calcControls['maxSSwave'][pfx]
[1618]	611	Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
[1630]	612	'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
[1618]	613	XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
	614	FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
	615	USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
	616	MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
[1630]	617	waveTypes = []
[1625]	618	keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
	619	'Tzero:':XSSdata[0],'Xslope:':XSSdata[1],'Yslope:':XSSdata[2],'Zslope:':XSSdata[3],
[1618]	620	'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
	621	'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
	622	'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
	623	'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
	624	for iatm in range(Natoms):
[1630]	625	waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
[1618]	626	for key in keys:
	627	for m in range(Nwave[key[0]]):
	628	parm = pfx+key+str(iatm)+':%d'%(m)
	629	if parm in parmDict:
[1630]	630	keys[key][m][iatm] = parmDict[parm]
	631	return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze()
[1618]	632
[1299]	633	def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
[1391]	634	''' Not Used: here only for comparison the StructureFactor2 - faster version
	635	Compute structure factors for all h,k,l for phase
[1299]	636	puts the result, F^2, in each ref[8] in refList
	637	input:
	638
	639	:param dict refDict: where
	640	'RefList' list where each ref = h,k,l,m,d,...
	641	'FF' dict of form factors - filed in below
	642	:param np.array G: reciprocal metric tensor
	643	:param str pfx: phase id string
	644	:param dict SGData: space group info. dictionary output from SpcGroup
	645	:param dict calcControls:
	646	:param dict ParmDict:
	647
	648	'''
	649	phfx = pfx.split(':')[0]+hfx
	650	ast = np.sqrt(np.diag(G))
	651	Mast = twopisq*np.multiply.outer(ast,ast)
	652	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
	653	SGT = np.array([ops[1] for ops in SGData['SGOps']])
	654	FFtables = calcControls['FFtables']
	655	BLtables = calcControls['BLtables']
	656	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
	657	FF = np.zeros(len(Tdata))
[1453]	658	if 'NC' in calcControls[hfx+'histType']:
[1391]	659	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
[1299]	660	else:
	661	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	662	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	663	Uij = np.array(G2lat.U6toUij(Uijdata))
	664	bij = Mast*Uij.T
	665	if not len(refDict['FF']):
	666	if 'N' in calcControls[hfx+'histType']:
[1360]	667	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
[1299]	668	else:
	669	dat = G2el.getFFvalues(FFtables,0.)
	670	refDict['FF']['El'] = dat.keys()
	671	refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
	672	for iref,refl in enumerate(refDict['RefList']):
[1453]	673	if 'NT' in calcControls[hfx+'histType']:
[1459]	674	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14])
[1299]	675	fbs = np.array([0,0])
	676	H = refl[:3]
	677	SQ = 1./(2.refl[4])*2
	678	SQfactor = 4.0SQtwopisq
	679	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
	680	if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor
	681	if 'N' in calcControls[hfx+'histType']:
	682	dat = G2el.getBLvalues(BLtables)
	683	refDict['FF']['FF'][iref] = dat.values()
	684	else: #'X'
	685	dat = G2el.getFFvalues(FFtables,SQ)
	686	refDict['FF']['FF'][iref] = dat.values()
	687	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
	688	FF = refDict['FF']['FF'][iref][Tindx]
	689	Uniq = np.inner(H,SGMT)
	690	Phi = np.inner(H,SGT)
	691	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
	692	sinp = np.sin(phase)
	693	cosp = np.cos(phase)
	694	biso = -SQfactor*Uisodata
	695	Tiso = np.where(biso<1.,np.exp(biso),1.0)
	696	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
	697	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
	698	Tcorr = TisoTuijMdata*Fdata/len(Uniq)
	699	fa = np.array([(FF+FP-Bab)cospTcorr,-FPPsinpTcorr])
	700	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
	701	if not SGData['SGInv']:
	702	fb = np.array([(FF+FP-Bab)sinpTcorr,FPPcospTcorr])
	703	fbs = np.sum(np.sum(fb,axis=1),axis=1)
	704	fasq = fas**2
	705	fbsq = fbs**2 #imaginary
	706	refl[9] = np.sum(fasq)+np.sum(fbsq)
	707	refl[10] = atan2d(fbs[0],fas[0])
	708
[1630]	709	def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
[1613]	710	'''
	711	Compute super structure factors for all h,k,l,m for phase
	712	puts the result, F^2, in each ref[8+im] in refList
	713	input:
	714
	715	:param dict refDict: where
	716	'RefList' list where each ref = h,k,l,m,d,...
	717	'FF' dict of form factors - filed in below
	718	:param np.array G: reciprocal metric tensor
	719	:param str pfx: phase id string
	720	:param dict SGData: space group info. dictionary output from SpcGroup
	721	:param dict calcControls:
	722	:param dict ParmDict:
	723
	724	'''
	725	phfx = pfx.split(':')[0]+hfx
	726	ast = np.sqrt(np.diag(G))
	727	Mast = twopisq*np.multiply.outer(ast,ast)
	728	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
	729	SGT = np.array([ops[1] for ops in SGData['SGOps']])
[1618]	730	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
	731	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
[1613]	732	FFtables = calcControls['FFtables']
	733	BLtables = calcControls['BLtables']
	734	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
[1630]	735	waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
[1613]	736	FF = np.zeros(len(Tdata))
	737	if 'NC' in calcControls[hfx+'histType']:
	738	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
	739	else:
	740	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	741	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	742	Uij = np.array(G2lat.U6toUij(Uijdata))
	743	bij = Mast*Uij.T
	744	if not len(refDict['FF']):
	745	if 'N' in calcControls[hfx+'histType']:
	746	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
	747	else:
	748	dat = G2el.getFFvalues(FFtables,0.)
	749	refDict['FF']['El'] = dat.keys()
	750	refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
	751	for iref,refl in enumerate(refDict['RefList']):
	752	if 'NT' in calcControls[hfx+'histType']:
	753	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
	754	fbs = np.array([0,0])
	755	H = refl[:4]
	756	SQ = 1./(2.refl[4+im])*2
	757	SQfactor = 4.0SQtwopisq
	758	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
	759	if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor
	760	if 'N' in calcControls[hfx+'histType']:
	761	dat = G2el.getBLvalues(BLtables)
	762	refDict['FF']['FF'][iref] = dat.values()
	763	else: #'X'
	764	dat = G2el.getFFvalues(FFtables,SQ)
	765	refDict['FF']['FF'][iref] = dat.values()
	766	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
	767	FF = refDict['FF']['FF'][iref][Tindx]
	768	Uniq = np.inner(H[:3],SGMT)
[1622]	769	SSUniq = np.inner(H,SSGMT)
[1613]	770	Phi = np.inner(H[:3],SGT)
[1622]	771	SSPhi = np.inner(H,SSGT)
[1630]	772	GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
[1613]	773	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
	774	sinp = np.sin(phase)
	775	cosp = np.cos(phase)
	776	biso = -SQfactor*Uisodata
	777	Tiso = np.where(biso<1.,np.exp(biso),1.0)
	778	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
	779	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
	780	Tcorr = TisoTuijMdata*Fdata/len(Uniq)
	781	fa = np.array([(FF+FP-Bab)cospTcorr,-FPPsinpTcorr])
[1635]	782	fb = np.zeros_like(fa)
[1613]	783	if not SGData['SGInv']:
	784	fb = np.array([(FF+FP-Bab)sinpTcorr,FPPcospTcorr])
[1635]	785	fa = faGfpuA-fbGfpuB
	786	fb = fbGfpuA+faGfpuB
	787	fas = np.real(np.sum(np.sum(fa,axis=1),axis=1)) #real
	788	fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
[1622]	789
[1613]	790	fasq = fas**2
	791	fbsq = fbs**2 #imaginary
	792	refl[9+im] = np.sum(fasq)+np.sum(fbsq)
	793	refl[10+im] = atan2d(fbs[0],fas[0])
	794
[1299]	795	def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
	796	''' Compute structure factors for all h,k,l for phase
	797	puts the result, F^2, in each ref[8] in refList
	798	input:
	799
	800	:param dict refDict: where
	801	'RefList' list where each ref = h,k,l,m,d,...
	802	'FF' dict of form factors - filed in below
	803	:param np.array G: reciprocal metric tensor
	804	:param str pfx: phase id string
	805	:param dict SGData: space group info. dictionary output from SpcGroup
	806	:param dict calcControls:
	807	:param dict ParmDict:
	808
	809	'''
	810	phfx = pfx.split(':')[0]+hfx
	811	ast = np.sqrt(np.diag(G))
	812	Mast = twopisq*np.multiply.outer(ast,ast)
	813	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
	814	SGT = np.array([ops[1] for ops in SGData['SGOps']])
	815	FFtables = calcControls['FFtables']
	816	BLtables = calcControls['BLtables']
[1878]	817	Flack = 1.0
	818	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']:
	819	Flack = 1.-2.*parmDict[phfx+'Flack']
[1886]	820	TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
	821	if 'S' in calcControls[hfx+'histType']:
	822	TwinLaw = calcControls[phfx+'TwinLaw']
	823	TwinFr = np.array([parmDict[phfx+'TwinFr;'+str(i)] for i in range(len(TwinLaw))])
	824	if len(TwinLaw) > 1:
	825	TwinFr[0] = 1.-np.sum(TwinFr[1:])
[1299]	826	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
	827	FF = np.zeros(len(Tdata))
[1453]	828	if 'NC' in calcControls[hfx+'histType']:
[1391]	829	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
[1453]	830	elif 'X' in calcControls[hfx+'histType']:
[1299]	831	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	832	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	833	Uij = np.array(G2lat.U6toUij(Uijdata))
	834	bij = Mast*Uij.T
	835	blkSize = 100 #no. of reflections in a block
	836	nRef = refDict['RefList'].shape[0]
	837	if not len(refDict['FF']): #no form factors - 1st time thru StructureFactor
	838	if 'N' in calcControls[hfx+'histType']:
	839	dat = G2el.getBLvalues(BLtables)
	840	refDict['FF']['El'] = dat.keys()
	841	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
	842	else: #'X'
	843	dat = G2el.getFFvalues(FFtables,0.)
	844	refDict['FF']['El'] = dat.keys()
	845	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
	846	for iref,ref in enumerate(refDict['RefList']):
	847	SQ = 1./(2.ref[4])*2
	848	dat = G2el.getFFvalues(FFtables,SQ)
	849	refDict['FF']['FF'][iref] *= dat.values()
	850	#reflection processing begins here - big arrays!
	851	iBeg = 0
	852	while iBeg < nRef:
	853	iFin = min(iBeg+blkSize,nRef)
[1886]	854	refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems)
	855	H = refl.T[:3] #array(blkSize,3)
	856	H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,3,nTwins) or (blkSize,3)
	857	SQ = 1./(2.refl.T[4])*2 #array(blkSize)
	858	SQfactor = 4.0SQtwopisq #ditto prev.
[1453]	859	if 'T' in calcControls[hfx+'histType']:
[1459]	860	if 'P' in calcControls[hfx+'histType']:
	861	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
	862	else:
	863	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
[1886]	864	FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
	865	FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
	866	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),len(SGT)*len(TwinLaw))
[1299]	867	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
[1886]	868	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
	869	Uniq = np.inner(H,SGMT)
	870	Phi = np.inner(H,SGT)
	871	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
[1299]	872	sinp = np.sin(phase)
	873	cosp = np.cos(phase)
	874	biso = -SQfactor*Uisodata[:,np.newaxis]
[1886]	875	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
	876	HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
[1299]	877	Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
[1886]	878	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/len(SGMT)
	879	if 'T' in calcControls[hfx+'histType']:
	880	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-np.reshape(FlackFPP,sinp.shape)sinp*Tcorr])
	881	fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinpTcorr,np.reshape(FlackFPP,cosp.shape)cosp*Tcorr])
	882	else:
	883	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-FlackFPPsinp*Tcorr])
	884	fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinpTcorr,FlackFPPcosp*Tcorr])
	885	fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl
	886	fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
[1786]	887	if SGData['SGInv']: #centrosymmetric; B=0
	888	fbs[0] *= 0.
	889	if 'P' in calcControls[hfx+'histType']:
	890	refl.T[9] = np.sum(fas2,axis=0)+np.sum(fbs2,axis=0)
[1886]	891	refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f"
[1786]	892	else:
[1886]	893	if len(TwinLaw) > 1:
	894	refl.T[9] = np.sum(fas[:,:,0]2,axis=0)+np.sum(fbs[:,:,0]2,axis=0) #FcT from primary twin element
	895	refl.T[7] = np.sum(TwinFrnp.sum(fas,axis=0)*2,axis=-1)+ \
	896	np.sum(TwinFrnp.sum(fbs,axis=0)*2,axis=-1) #Fc sum over twins
	897	# what goes in refl.T[8]? (FoT)
	898	refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f"
	899	else:
	900	refl.T[9] = np.sum(fas,axis=0)2+np.sum(fbs,axis=0)2
	901	refl.T[7] = np.copy(refl.T[9])
	902	refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f"
[1786]	903	# refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
[1299]	904	iBeg += blkSize
	905
	906	def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
	907	'Needs a doc string'
	908	phfx = pfx.split(':')[0]+hfx
	909	ast = np.sqrt(np.diag(G))
	910	Mast = twopisq*np.multiply.outer(ast,ast)
	911	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
	912	SGT = np.array([ops[1] for ops in SGData['SGOps']])
	913	FFtables = calcControls['FFtables']
	914	BLtables = calcControls['BLtables']
[1886]	915	TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
	916	if 'S' in calcControls[hfx+'histType']:
	917	TwinLaw = calcControls[phfx+'TwinLaw']
	918	TwinFr = np.array([parmDict[phfx+'TwinFr;'+str(i)] for i in range(len(TwinLaw))])
	919	if len(TwinLaw) > 1:
[1887]	920	TwinFr[0] = 1.-np.sum(TwinFr[1:])
	921	nTwin = len(TwinLaw)
[1299]	922	nRef = len(refDict['RefList'])
	923	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
	924	mSize = len(Mdata)
	925	FF = np.zeros(len(Tdata))
[1459]	926	if 'NC' in calcControls[hfx+'histType']:
[1391]	927	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
[1459]	928	elif 'X' in calcControls[hfx+'histType']:
[1299]	929	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	930	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	931	Uij = np.array(G2lat.U6toUij(Uijdata))
	932	bij = Mast*Uij.T
	933	dFdvDict = {}
[1887]	934	dFdfr = np.squeeze(np.zeros((nRef,nTwin,mSize)))
	935	dFdx = np.squeeze(np.zeros((nRef,nTwin,mSize,3)))
	936	dFdui = np.squeeze(np.zeros((nRef,nTwin,mSize)))
	937	dFdua = np.squeeze(np.zeros((nRef,nTwin,mSize,6)))
	938	dFdbab = np.squeeze(np.zeros((nRef,nTwin,2)))
	939	dFdfl = np.squeeze(np.zeros((nRef,nTwin)))
	940	dFdtw = np.zeros((nRef,nTwin))
[1878]	941	Flack = 1.0
	942	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']:
	943	Flack = 1.-2.*parmDict[phfx+'Flack']
[1299]	944	for iref,refl in enumerate(refDict['RefList']):
[1459]	945	if 'T' in calcControls[hfx+'histType']:
	946	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
[1299]	947	H = np.array(refl[:3])
[1888]	948	H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)
[1299]	949	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
	950	SQfactor = 8.0SQnp.pi**2
	951	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
	952	Bab = parmDict[phfx+'BabA']*dBabdA
	953	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
[1886]	954	FF = refDict['FF']['FF'][iref].T[Tindx].T
[1888]	955	Uniq = np.inner(H,SGMT) # array(3,nSGOp) or
[1299]	956	Phi = np.inner(H,SGT)
[1886]	957	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
[1299]	958	sinp = np.sin(phase)
	959	cosp = np.cos(phase)
	960	occ = Mdata*Fdata/len(Uniq)
[1886]	961	biso = -SQfactor*Uisodata[:,np.newaxis]
[1887]	962	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1).T
[1886]	963	HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
	964	Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
[1887]	965	Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
[1886]	966	Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
	967	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/len(SGMT)
[1299]	968	fot = (FF+FP-Bab)occTcorr
[1886]	969	fotp = FPPoccTcorr
[1887]	970	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FlackFPPsinp*Tcorr])
	971	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FlackFPPcosp*Tcorr])
[1886]	972	fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl
	973	fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
	974	fax = np.array([-fotsinp,-fotpcosp]) #positions
	975	fbx = np.array([fotcosp,-fotpsinp])
[1887]	976	#sum below is over Uniq
[1886]	977	dfadfr = np.sum(fa/occ,axis=-2) #Fdata != 0 ever avoids /0. problem
[1887]	978	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
[1886]	979	dfadui = np.sum(-SQfactor*fa,axis=-2)
[1887]	980	if len(TwinLaw) > 1:
[1888]	981	dfadx = np.array([np.sum(twopiUniq[it]np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
	982	dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
[1887]	983	else:
	984	dfadx = np.sum(twopiUniqnp.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
	985	dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
[1299]	986	if not SGData['SGInv']:
[1886]	987	dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 ever avoids /0. problem
[1299]	988	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
[1879]	989	dfadfl = np.sum(-fotp[:,np.newaxis]*sinp)
	990	dfbdfl = np.sum(fotp[:,np.newaxis]*cosp)
[1887]	991	dfbdui = np.sum(-SQfactor*fb,axis=-2)
	992	if len(TwinLaw) > 1:
[1888]	993	dfbdx = np.array([np.sum(twopiUniq[it]np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
	994	dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
[1887]	995	else:
	996	dfbdx = np.sum(twopiUniqnp.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)
	997	dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
[1864]	998	else:
	999	dfbdfr = np.zeros_like(dfadfr)
	1000	dfbdx = np.zeros_like(dfadx)
	1001	dfbdui = np.zeros_like(dfadui)
	1002	dfbdua = np.zeros_like(dfadua)
	1003	dfbdba = np.zeros_like(dfadba)
[1878]	1004	dfadfl = 0.0
	1005	dfbdfl = 0.0
[1864]	1006	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!
	1007	if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
	1008	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)+ \
	1009	2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
	1010	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])+ \
	1011	2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
	1012	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])+ \
	1013	2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
	1014	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])+ \
	1015	2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
	1016	else:
	1017	SA = fas[0]-fbs[1]
	1018	SB = fbs[0]+fas[1]
[1887]	1019	if nTwin > 1:
[1888]	1020	dFdfr[iref] = [2.SA[it](dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \
	1021	2.SB[it](dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
	1022	dFdx[iref] = [2.SA[it](dfadx[0][it]+dfbdx[1][it])+2.SB[it](dfbdx[0][it]+dfadx[1][it]) for it in range(nTwin)]
	1023	dFdui[iref] = [2.SA[it](dfadui[0][it]+dfbdui[1][it])+2.SB[it](dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]
[1887]	1024	dFdua[iref] = [2.SA[it](dfadua[it][0]+dfbdua[it][1])+2.SB[it](dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
	1025	dFdfl[iref] = -SAdfadfl-SBdfbdfl
	1026	dFdtw[iref] = 2.SA+2.SB
	1027	else:
	1028	dFdfr[iref] = 2.SA(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
	1029	2.SB(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq)
	1030	dFdx[iref] = 2.SA(dfadx[0]+dfbdx[1])+2.SB(dfbdx[0]+dfadx[1])
	1031	dFdui[iref] = 2.SA(dfadui[0]+dfbdui[1])+2.SB(dfbdui[0]+dfadui[1])
	1032	dFdua[iref] = 2.SA(dfadua[0]+dfbdua[1])+2.SB(dfbdua[0]+dfadua[1])
	1033	dFdfl[iref] = -SAdfadfl-SBdfbdfl
[1864]	1034	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T+ \
	1035	2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
[1878]	1036
[1299]	1037	#loop over atoms - each dict entry is list of derivatives for all the reflections
[1887]	1038	if nTwin > 1:
	1039	for i in range(len(Mdata)):
	1040	dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
	1041	dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
	1042	dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
	1043	dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
	1044	dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
	1045	dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
	1046	dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
	1047	dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
	1048	dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5dFdua.T[3][i]TwinFr[:,np.newaxis],axis=0)
	1049	dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5dFdua.T[4][i]TwinFr[:,np.newaxis],axis=0)
	1050	dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5dFdua.T[5][i]TwinFr[:,np.newaxis],axis=0)
	1051	dFdvDict[phfx+'Flack'] = np.sum(dFdfl.T*TwinFr[:,np.newaxis],axis=0)
	1052	else:
	1053	for i in range(len(Mdata)):
	1054	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
	1055	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
	1056	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
	1057	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
	1058	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
	1059	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
	1060	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
	1061	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
	1062	dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
	1063	dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
	1064	dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
	1065	dFdvDict[phfx+'Flack'] = dFdfl.T
[1878]	1066	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
	1067	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
[1887]	1068	if nTwin > 1:
	1069	for i in range(nTwin-1): #skip the base twin element
	1070	dFdvDict[phfx+'TwinFr;'+str(i+1)] = dFdtw.T[i+1]
[1299]	1071	return dFdvDict
	1072
[1630]	1073	def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
[1613]	1074	'Needs a doc string'
	1075	phfx = pfx.split(':')[0]+hfx
	1076	ast = np.sqrt(np.diag(G))
	1077	Mast = twopisq*np.multiply.outer(ast,ast)
	1078	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
	1079	SGT = np.array([ops[1] for ops in SGData['SGOps']])
[1630]	1080	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
	1081	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
[1613]	1082	FFtables = calcControls['FFtables']
	1083	BLtables = calcControls['BLtables']
	1084	nRef = len(refDict['RefList'])
	1085	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
[1630]	1086	waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
[1613]	1087	mSize = len(Mdata)
	1088	FF = np.zeros(len(Tdata))
	1089	if 'NC' in calcControls[hfx+'histType']:
	1090	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
	1091	elif 'X' in calcControls[hfx+'histType']:
	1092	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	1093	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	1094	Uij = np.array(G2lat.U6toUij(Uijdata))
	1095	bij = Mast*Uij.T
	1096	dFdvDict = {}
	1097	dFdfr = np.zeros((nRef,mSize))
	1098	dFdx = np.zeros((nRef,mSize,3))
	1099	dFdui = np.zeros((nRef,mSize))
	1100	dFdua = np.zeros((nRef,mSize,6))
	1101	dFdbab = np.zeros((nRef,2))
	1102	for iref,refl in enumerate(refDict['RefList']):
	1103	if 'T' in calcControls[hfx+'histType']:
	1104	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
	1105	H = np.array(refl[:4])
	1106	SQ = 1./(2.refl[4+im])2 # or (sin(theta)/lambda)*2
	1107	SQfactor = 8.0SQnp.pi**2
	1108	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
	1109	Bab = parmDict[phfx+'BabA']*dBabdA
	1110	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
	1111	FF = refDict['FF']['FF'][iref].T[Tindx]
[1630]	1112	Uniq = np.inner(H[:3],SGMT)
	1113	SSUniq = np.inner(H,SSGMT)
	1114	Phi = np.inner(H[:3],SGT)
	1115	SSPhi = np.inner(H,SSGT)
[1635]	1116	GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
	1117	dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
	1118	#need ModulationDerv here dGAdXsin, etc
[1613]	1119	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
	1120	sinp = np.sin(phase)
	1121	cosp = np.cos(phase)
	1122	occ = Mdata*Fdata/len(Uniq)
	1123	biso = -SQfactor*Uisodata
	1124	Tiso = np.where(biso<1.,np.exp(biso),1.0)
[1635]	1125	HbH = -np.inner(H[:3],np.inner(bij,H[:3]))
[1613]	1126	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
	1127	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
	1128	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
	1129	Tcorr = Tiso*Tuij
	1130	fot = (FF+FP-Bab)occTcorr
	1131	fotp = FPPoccTcorr
	1132	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
	1133	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
[1635]	1134	GfpuA = np.swapaxes(GfpuA,1,2)
	1135	GfpuB = np.swapaxes(GfpuB,1,2)
	1136	fa = faGfpuA-fbGfpuB
	1137	fb = fbGfpuA+faGfpuB
[1613]	1138
	1139	fas = np.sum(np.sum(fa,axis=1),axis=1)
	1140	fbs = np.sum(np.sum(fb,axis=1),axis=1)
	1141	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
	1142	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
[1635]	1143	fax = faxGfpuA-fbxGfpuB
	1144	fbx = fbxGfpuA+faxGfpuB
[1613]	1145	#sum below is over Uniq
	1146	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem
	1147	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
	1148	dfadui = np.sum(-SQfactor*fa,axis=2)
	1149	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
	1150	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
	1151	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!
	1152	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
	1153	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
	1154	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
	1155	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
	1156	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
[1635]	1157	#need dFdXsin, etc for modulations...
[1613]	1158	if not SGData['SGInv']:
	1159	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem
	1160	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
	1161	dfbdui = np.sum(-SQfactor*fb,axis=2)
	1162	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
	1163	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
	1164	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
	1165	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
	1166	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
	1167	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
	1168	dFdbab[iref] += 2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
[1635]	1169	#need dFdXsin, etc for modulations...
[1613]	1170	#loop over atoms - each dict entry is list of derivatives for all the reflections
	1171	for i in range(len(Mdata)):
	1172	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
	1173	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
	1174	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
	1175	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
	1176	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
	1177	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
	1178	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
	1179	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
[1857]	1180	dFdvDict[pfx+'AU12:'+str(i)] = .5*dFdua.T[3][i]
	1181	dFdvDict[pfx+'AU13:'+str(i)] = .5*dFdua.T[4][i]
	1182	dFdvDict[pfx+'AU23:'+str(i)] = .5*dFdua.T[5][i]
[1635]	1183	#need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
[1878]	1184	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
	1185	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
[1613]	1186	return dFdvDict
	1187
[1595]	1188	def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
[1465]	1189	''' Single crystal extinction function; returns extinction & derivative
[1299]	1190	'''
[1460]	1191	extCor = 1.0
[1299]	1192	dervDict = {}
[1857]	1193	dervCor = 1.0
[1299]	1194	if calcControls[phfx+'EType'] != 'None':
[1595]	1195	SQ = 1/(4.ref[4+im]*2)
[1453]	1196	if 'C' in parmDict[hfx+'Type']:
	1197	cos2T = 1.0-2.SQparmDict[hfx+'Lam']**2 #cos(2theta)
	1198	else: #'T'
[1595]	1199	cos2T = 1.0-2.SQref[12+im]**2 #cos(2theta)
[1299]	1200	if 'SXC' in parmDict[hfx+'Type']:
	1201	AV = 7.9406e5/parmDict[pfx+'Vol']**2
	1202	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
	1203	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
[1776]	1204	PLZ = AVP12ref[9+im]parmDict[hfx+'Lam']*2
[1299]	1205	elif 'SNT' in parmDict[hfx+'Type']:
	1206	AV = 1.e7/parmDict[pfx+'Vol']**2
[1453]	1207	PL = SQ
[1776]	1208	PLZ = AVref[9+im]ref[12+im]**2
[1299]	1209	elif 'SNC' in parmDict[hfx+'Type']:
	1210	AV = 1.e7/parmDict[pfx+'Vol']**2
	1211	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
[1776]	1212	PLZ = AVref[9+im]parmDict[hfx+'Lam']**2
[1299]	1213
	1214	if 'Primary' in calcControls[phfx+'EType']:
	1215	PLZ *= 1.5
	1216	else:
[1453]	1217	if 'C' in parmDict[hfx+'Type']:
	1218	PLZ *= calcControls[phfx+'Tbar']
	1219	else: #'T'
[1595]	1220	PLZ *= ref[13+im] #t-bar
[1299]	1221	if 'Primary' in calcControls[phfx+'EType']:
[1453]	1222	PLZ *= 1.5
[1299]	1223	PSIG = parmDict[phfx+'Ep']
	1224	elif 'I & II' in calcControls[phfx+'EType']:
	1225	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
	1226	elif 'Type II' in calcControls[phfx+'EType']:
	1227	PSIG = parmDict[phfx+'Es']
	1228	else: # 'Secondary Type I'
	1229	PSIG = parmDict[phfx+'Eg']/PL
	1230
	1231	AG = 0.58+0.48cos2T+0.24cos2T**2
	1232	AL = 0.025+0.285*cos2T
	1233	BG = 0.02-0.025*cos2T
	1234	BL = 0.15-0.2(0.75-cos2T)*2
	1235	if cos2T < 0.:
	1236	BL = -0.45*cos2T
	1237	CG = 2.
	1238	CL = 2.
	1239	PF = PLZ*PSIG
	1240
	1241	if 'Gaussian' in calcControls[phfx+'EApprox']:
	1242	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
	1243	extCor = np.sqrt(PF4)
	1244	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
	1245	else:
	1246	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
	1247	extCor = np.sqrt(PF4)
	1248	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
	1249
[1857]	1250	dervCor = (1.+PF)*PF3 #extinction corr for other derivatives
[1299]	1251	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
[1595]	1252	dervDict[phfx+'Ep'] = -ref[7+im]PLZPF3
[1299]	1253	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
[1595]	1254	dervDict[phfx+'Es'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
[1299]	1255	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
[1595]	1256	dervDict[phfx+'Eg'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
[1299]	1257
[1857]	1258	return 1./extCor,dervDict,dervCor
[1299]	1259
	1260	def Dict2Values(parmdict, varylist):
	1261	'''Use before call to leastsq to setup list of values for the parameters
	1262	in parmdict, as selected by key in varylist'''
	1263	return [parmdict[key] for key in varylist]
	1264
[1320]	1265	def Values2Dict(parmdict, varylist, values):
	1266	''' Use after call to leastsq to update the parameter dictionary with
	1267	values corresponding to keys in varylist'''
	1268	parmdict.update(zip(varylist,values))
[1299]	1269
	1270	def GetNewCellParms(parmDict,varyList):
	1271	'Needs a doc string'
	1272	newCellDict = {}
	1273	Anames = ['A'+str(i) for i in range(6)]
	1274	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
	1275	for item in varyList:
	1276	keys = item.split(':')
	1277	if keys[2] in Ddict:
	1278	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
	1279	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
[1562]	1280	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
[1482]	1281	return newCellDict # is e.g. {'0::D11':A0-D11}
[1299]	1282
	1283	def ApplyXYZshifts(parmDict,varyList):
	1284	'''
	1285	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
	1286
	1287	:param dict parmDict: parameter dictionary
[1335]	1288	:param list varyList: list of variables (not used!)
[1299]	1289	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
	1290
	1291	'''
	1292	newAtomDict = {}
	1293	for item in parmDict:
	1294	if 'dA' in item:
	1295	parm = ''.join(item.split('d'))
	1296	parmDict[parm] += parmDict[item]
	1297	newAtomDict[item] = [parm,parmDict[parm]]
	1298	return newAtomDict
	1299
[1595]	1300	def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
[1299]	1301	'Spherical harmonics texture'
	1302	IFCoup = 'Bragg' in calcControls[hfx+'instType']
[1459]	1303	if 'T' in calcControls[hfx+'histType']:
	1304	tth = parmDict[hfx+'2-theta']
	1305	else:
[1595]	1306	tth = refl[5+im]
[1299]	1307	odfCor = 1.0
	1308	H = refl[:3]
	1309	cell = G2lat.Gmat2cell(g)
	1310	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
[1767]	1311	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
[1299]	1312	phi,beta = G2lat.CrsAng(H,cell,SGData)
[1459]	1313	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
[1299]	1314	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
	1315	for item in SHnames:
	1316	L,M,N = eval(item.strip('C'))
	1317	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
	1318	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
	1319	Lnorm = G2lat.Lnorm(L)
	1320	odfCor += parmDict[pfx+item]LnormKcl*Ksl
	1321	return odfCor
	1322
[1595]	1323	def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
[1299]	1324	'Spherical harmonics texture derivatives'
[1459]	1325	if 'T' in calcControls[hfx+'histType']:
	1326	tth = parmDict[hfx+'2-theta']
	1327	else:
[1595]	1328	tth = refl[5+im]
[1299]	1329	IFCoup = 'Bragg' in calcControls[hfx+'instType']
	1330	odfCor = 1.0
	1331	dFdODF = {}
	1332	dFdSA = [0,0,0]
	1333	H = refl[:3]
	1334	cell = G2lat.Gmat2cell(g)
	1335	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
[1767]	1336	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
[1299]	1337	phi,beta = G2lat.CrsAng(H,cell,SGData)
[1459]	1338	psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
[1299]	1339	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
	1340	for item in SHnames:
	1341	L,M,N = eval(item.strip('C'))
	1342	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
	1343	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
	1344	Lnorm = G2lat.Lnorm(L)
	1345	odfCor += parmDict[pfx+item]LnormKcl*Ksl
	1346	dFdODF[pfx+item] = LnormKclKsl
	1347	for i in range(3):
	1348	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
	1349	return odfCor,dFdODF,dFdSA
	1350
[1595]	1351	def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
[1299]	1352	'spherical harmonics preferred orientation (cylindrical symmetry only)'
[1459]	1353	if 'T' in calcControls[hfx+'histType']:
	1354	tth = parmDict[hfx+'2-theta']
	1355	else:
[1595]	1356	tth = refl[5+im]
[1299]	1357	odfCor = 1.0
	1358	H = refl[:3]
	1359	cell = G2lat.Gmat2cell(g)
[1835]	1360	Sangls = [0.,0.,0.]
[1299]	1361	if 'Bragg' in calcControls[hfx+'instType']:
	1362	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
	1363	IFCoup = True
	1364	else:
[1791]	1365	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
[1299]	1366	IFCoup = False
	1367	phi,beta = G2lat.CrsAng(H,cell,SGData)
[1835]	1368	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
[1777]	1369	SHnames = calcControls[phfx+'SHnames']
[1299]	1370	for item in SHnames:
	1371	L,N = eval(item.strip('C'))
[1812]	1372	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
	1373	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
	1374	Lnorm = G2lat.Lnorm(L)
	1375	odfCor += parmDict[phfx+item]LnormKcl*Ksl
[1299]	1376	return np.squeeze(odfCor)
	1377
[1595]	1378	def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
[1299]	1379	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
[1459]	1380	if 'T' in calcControls[hfx+'histType']:
	1381	tth = parmDict[hfx+'2-theta']
	1382	else:
[1595]	1383	tth = refl[5+im]
[1299]	1384	odfCor = 1.0
	1385	dFdODF = {}
	1386	H = refl[:3]
	1387	cell = G2lat.Gmat2cell(g)
[1835]	1388	Sangls = [0.,0.,0.]
[1299]	1389	if 'Bragg' in calcControls[hfx+'instType']:
	1390	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
	1391	IFCoup = True
	1392	else:
[1791]	1393	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
[1299]	1394	IFCoup = False
	1395	phi,beta = G2lat.CrsAng(H,cell,SGData)
[1835]	1396	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
[1777]	1397	SHnames = calcControls[phfx+'SHnames']
[1299]	1398	for item in SHnames:
	1399	L,N = eval(item.strip('C'))
[1812]	1400	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
	1401	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
	1402	Lnorm = G2lat.Lnorm(L)
	1403	odfCor += parmDict[phfx+item]LnormKcl*Ksl
	1404	dFdODF[phfx+item] = KclKslLnorm
[1299]	1405	return odfCor,dFdODF
	1406
[1595]	1407	def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
[1460]	1408	'March-Dollase preferred orientation correction'
[1299]	1409	POcorr = 1.0
[1460]	1410	MD = parmDict[phfx+'MD']
	1411	if MD != 1.0:
	1412	MDAxis = calcControls[phfx+'MDAxis']
	1413	sumMD = 0
	1414	for H in uniq:
	1415	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
	1416	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
	1417	sumMD += A**3
	1418	POcorr = sumMD/len(uniq)
[1299]	1419	return POcorr
	1420
[1595]	1421	def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
[1299]	1422	'Needs a doc string'
	1423	POcorr = 1.0
	1424	POderv = {}
	1425	if calcControls[phfx+'poType'] == 'MD':
	1426	MD = parmDict[phfx+'MD']
	1427	MDAxis = calcControls[phfx+'MDAxis']
	1428	sumMD = 0
	1429	sumdMD = 0
	1430	for H in uniq:
	1431	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
	1432	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
	1433	sumMD += A**3
	1434	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
	1435	POcorr = sumMD/len(uniq)
	1436	POderv[phfx+'MD'] = sumdMD/len(uniq)
	1437	else: #spherical harmonics
	1438	if calcControls[phfx+'SHord']:
[1595]	1439	POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
[1299]	1440	return POcorr,POderv
	1441
[1595]	1442	def GetAbsorb(refl,im,hfx,calcControls,parmDict):
[1299]	1443	'Needs a doc string'
	1444	if 'Debye' in calcControls[hfx+'instType']:
[1459]	1445	if 'T' in calcControls[hfx+'histType']:
[1766]	1446	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
[1459]	1447	else:
[1595]	1448	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
[1299]	1449	else:
[1595]	1450	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
[1299]	1451
[1595]	1452	def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
[1299]	1453	'Needs a doc string'
	1454	if 'Debye' in calcControls[hfx+'instType']:
[1459]	1455	if 'T' in calcControls[hfx+'histType']:
[1766]	1456	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
[1459]	1457	else:
[1595]	1458	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
[1299]	1459	else:
[1595]	1460	return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
[1460]	1461
[1595]	1462	def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
[1460]	1463	'Needs a doc string'
	1464	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
[1465]	1465	pi2 = np.sqrt(2./np.pi)
[1460]	1466	if 'T' in calcControls[hfx+'histType']:
[1766]	1467	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
[1595]	1468	wave = refl[14+im]
[1460]	1469	else: #'C'W
[1595]	1470	sth2 = sind(refl[5+im]/2.)**2
[1460]	1471	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
	1472	c2th = 1.-2.0*sth2
[1595]	1473	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
[1460]	1474	if 'X' in calcControls[hfx+'histType']:
	1475	flv2 = 0.079411(1.0+c2th**2)/2.0
	1476	xfac = flv2*parmDict[phfx+'Extinction']
	1477	exb = 1.0
	1478	if xfac > -1.:
[1874]	1479	exb = 1./np.sqrt(1.+xfac)
[1460]	1480	exl = 1.0
	1481	if 0 < xfac <= 1.:
	1482	xn = np.array([xfac**(i+1) for i in range(6)])
[1761]	1483	exl += np.sum(xn*coef)
[1460]	1484	elif xfac > 1.:
	1485	xfac2 = 1./np.sqrt(xfac)
	1486	exl = pi2(1.-0.125/xfac)xfac2
	1487	return exbsth2+exl(1.-sth2)
[1299]	1488
[1595]	1489	def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
[1460]	1490	'Needs a doc string'
[1761]	1491	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
	1492	pi2 = np.sqrt(2./np.pi)
	1493	if 'T' in calcControls[hfx+'histType']:
[1766]	1494	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
[1761]	1495	wave = refl[14+im]
	1496	else: #'C'W
	1497	sth2 = sind(refl[5+im]/2.)**2
	1498	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
	1499	c2th = 1.-2.0*sth2
	1500	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
	1501	if 'X' in calcControls[hfx+'histType']:
	1502	flv2 = 0.079411(1.0+c2th**2)/2.0
	1503	xfac = flv2*parmDict[phfx+'Extinction']
	1504	dbde = -500.*flv2
	1505	if xfac > -1.:
[1874]	1506	dbde = -0.5flv2/np.sqrt(1.+xfac)*3
[1761]	1507	dlde = 0.
	1508	if 0 < xfac <= 1.:
	1509	xn = np.array([iflv2xfac**i for i in [1,2,3,4,5,6]])
	1510	dlde = np.sum(xn*coef)
	1511	elif xfac > 1.:
	1512	xfac2 = 1./np.sqrt(xfac)
	1513	dlde = flv2pi2xfac2(-1./xfac+0.375/xfac*2)
	1514
	1515	return dbdesth2+dlde(1.-sth2)
[1460]	1516
[1595]	1517	def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
[1299]	1518	'Needs a doc string' #need powder extinction!
[1595]	1519	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3+im] #scale*multiplicity
[1299]	1520	if 'X' in parmDict[hfx+'Type']:
[1595]	1521	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
[1460]	1522	POcorr = 1.0
[1770]	1523	if pfx+'SHorder' in parmDict: #generalized spherical harmonics texture - takes precidence
[1595]	1524	POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
[1460]	1525	elif calcControls[phfx+'poType'] == 'MD': #March-Dollase
[1595]	1526	POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
[1460]	1527	elif calcControls[phfx+'SHord']: #cylindrical spherical harmonics
[1595]	1528	POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
[1460]	1529	Icorr *= POcorr
	1530	AbsCorr = 1.0
[1595]	1531	AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
[1460]	1532	Icorr *= AbsCorr
[1595]	1533	ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
[1460]	1534	Icorr *= ExtCorr
	1535	return Icorr,POcorr,AbsCorr,ExtCorr
[1299]	1536
[1595]	1537	def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
[1299]	1538	'Needs a doc string' #need powder extinction derivs!
	1539	dIdsh = 1./parmDict[hfx+'Scale']
	1540	dIdsp = 1./parmDict[phfx+'Scale']
	1541	if 'X' in parmDict[hfx+'Type']:
[1595]	1542	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
[1299]	1543	dIdPola /= pola
	1544	else: #'N'
	1545	dIdPola = 0.0
	1546	dFdODF = {}
	1547	dFdSA = [0,0,0]
[1461]	1548	dIdPO = {}
[1299]	1549	if pfx+'SHorder' in parmDict:
[1595]	1550	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
[1299]	1551	for iSH in dFdODF:
	1552	dFdODF[iSH] /= odfCor
	1553	for i in range(3):
	1554	dFdSA[i] /= odfCor
[1460]	1555	elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
[1595]	1556	POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
[1460]	1557	for iPO in dIdPO:
	1558	dIdPO[iPO] /= POcorr
[1466]	1559	if 'T' in parmDict[hfx+'Type']:
[1595]	1560	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
	1561	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im] #/ext corr
[1466]	1562	else:
[1595]	1563	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
	1564	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im] #/ext corr
[1460]	1565	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
[1299]	1566
[1595]	1567	def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
[1299]	1568	'Needs a doc string'
[1484]	1569	if 'C' in calcControls[hfx+'histType']: #All checked & OK
[1595]	1570	costh = cosd(refl[5+im]/2.)
[1459]	1571	#crystallite size
	1572	if calcControls[phfx+'SizeType'] == 'isotropic':
	1573	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
	1574	elif calcControls[phfx+'SizeType'] == 'uniaxial':
	1575	H = np.array(refl[:3])
	1576	P = np.array(calcControls[phfx+'SizeAxis'])
	1577	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1578	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
	1579	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
	1580	else: #ellipsoidal crystallites
	1581	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
	1582	H = np.array(refl[:3])
	1583	lenR = G2pwd.ellipseSize(H,Sij,GB)
	1584	Sgam = 1.8wave/(np.picosth*lenR)
	1585	#microstrain
	1586	if calcControls[phfx+'MustrainType'] == 'isotropic':
[1595]	1587	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
[1459]	1588	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
	1589	H = np.array(refl[:3])
	1590	P = np.array(calcControls[phfx+'MustrainAxis'])
	1591	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1592	Si = parmDict[phfx+'Mustrain;i']
	1593	Sa = parmDict[phfx+'Mustrain;a']
[1595]	1594	Mgam = 0.018SiSatand(refl[5+im]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
[1459]	1595	else: #generalized - P.W. Stephens model
[1479]	1596	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
	1597	Sum = 0
	1598	for i,strm in enumerate(Strms):
	1599	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
[1595]	1600	Mgam = 0.018refl[4+im]2tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
[1484]	1601	elif 'T' in calcControls[hfx+'histType']: #All checked & OK
[1459]	1602	#crystallite size
[1484]	1603	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
[1595]	1604	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
[1484]	1605	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
[1459]	1606	H = np.array(refl[:3])
	1607	P = np.array(calcControls[phfx+'SizeAxis'])
	1608	cosP,sinP = G2lat.CosSinAngle(H,P,G)
[1595]	1609	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]*2/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a'])
[1484]	1610	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
	1611	else: #ellipsoidal crystallites #OK
[1459]	1612	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
	1613	H = np.array(refl[:3])
	1614	lenR = G2pwd.ellipseSize(H,Sij,GB)
[1595]	1615	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/lenR
[1459]	1616	#microstrain
[1484]	1617	if calcControls[phfx+'MustrainType'] == 'isotropic': #OK
[1595]	1618	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
[1484]	1619	elif calcControls[phfx+'MustrainType'] == 'uniaxial': #OK
[1459]	1620	H = np.array(refl[:3])
	1621	P = np.array(calcControls[phfx+'MustrainAxis'])
	1622	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1623	Si = parmDict[phfx+'Mustrain;i']
	1624	Sa = parmDict[phfx+'Mustrain;a']
[1595]	1625	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa/np.sqrt((SicosP)2+(SasinP)**2)
[1484]	1626	else: #generalized - P.W. Stephens model OK
	1627	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
[1479]	1628	Sum = 0
	1629	for i,strm in enumerate(Strms):
	1630	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
[1595]	1631	Mgam = 1.e-6parmDict[hfx+'difC']np.sqrt(Sum)refl[4+im]*3
[1459]	1632
[1299]	1633	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
	1634	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
	1635	sig /= ateln2
	1636	return sig,gam
	1637
[1595]	1638	def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
[1299]	1639	'Needs a doc string'
	1640	gamDict = {}
	1641	sigDict = {}
[1484]	1642	if 'C' in calcControls[hfx+'histType']: #All checked & OK
[1595]	1643	costh = cosd(refl[5+im]/2.)
	1644	tanth = tand(refl[5+im]/2.)
[1459]	1645	#crystallite size derivatives
	1646	if calcControls[phfx+'SizeType'] == 'isotropic':
	1647	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
	1648	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
	1649	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
	1650	elif calcControls[phfx+'SizeType'] == 'uniaxial':
	1651	H = np.array(refl[:3])
	1652	P = np.array(calcControls[phfx+'SizeAxis'])
	1653	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1654	Si = parmDict[phfx+'Size;i']
	1655	Sa = parmDict[phfx+'Size;a']
[1484]	1656	gami = 1.8wave/(costhnp.piSiSa)
[1459]	1657	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
	1658	Sgam = gami*sqtrm
[1484]	1659	dsi = gamiSicosP**2/sqtrm-Sgam/Si
	1660	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
[1459]	1661	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
	1662	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
	1663	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1664	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1665	else: #ellipsoidal crystallites
	1666	const = 1.8wave/(np.picosth)
	1667	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
	1668	H = np.array(refl[:3])
	1669	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
	1670	Sgam = const/lenR
	1671	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
	1672	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
	1673	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1674	gamDict[phfx+'Size;mx'] = Sgam
	1675	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
	1676
	1677	#microstrain derivatives
	1678	if calcControls[phfx+'MustrainType'] == 'isotropic':
[1595]	1679	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
[1459]	1680	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
	1681	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
	1682	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
	1683	H = np.array(refl[:3])
	1684	P = np.array(calcControls[phfx+'MustrainAxis'])
	1685	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1686	Si = parmDict[phfx+'Mustrain;i']
	1687	Sa = parmDict[phfx+'Mustrain;a']
	1688	gami = 0.018SiSa*tanth/np.pi
	1689	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
	1690	Mgam = gami/sqtrm
	1691	dsi = -gamiSicosP2/sqtrm3
	1692	dsa = -gamiSasinP2/sqtrm3
	1693	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
	1694	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
	1695	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1696	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1697	else: #generalized - P.W. Stephens model
[1595]	1698	const = 0.018refl[4+im]2tanth/np.pi
[1479]	1699	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
	1700	Sum = 0
	1701	for i,strm in enumerate(Strms):
	1702	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
	1703	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
	1704	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
	1705	Mgam = const*np.sqrt(Sum)
	1706	for i in range(len(Strms)):
	1707	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
	1708	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
[1459]	1709	gamDict[phfx+'Mustrain;mx'] = Mgam
	1710	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
[1484]	1711	else: #'T'OF - All checked & OK
	1712	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
[1595]	1713	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
[1484]	1714	gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
	1715	sigDict[phfx+'Size;i'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])*2/(ateln2parmDict[phfx+'Size;i'])
	1716	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
[1595]	1717	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
[1459]	1718	H = np.array(refl[:3])
	1719	P = np.array(calcControls[phfx+'SizeAxis'])
	1720	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1721	Si = parmDict[phfx+'Size;i']
	1722	Sa = parmDict[phfx+'Size;a']
[1484]	1723	gami = const/(Si*Sa)
[1459]	1724	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
[1484]	1725	Sgam = gami*sqtrm
	1726	dsi = gamiSicosP**2/sqtrm-Sgam/Si
	1727	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
	1728	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
	1729	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
[1459]	1730	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1731	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
[1484]	1732	else: #OK ellipsoidal crystallites
[1595]	1733	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
[1459]	1734	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
	1735	H = np.array(refl[:3])
	1736	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
	1737	Sgam = const/lenR
	1738	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
	1739	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
	1740	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
[1484]	1741	gamDict[phfx+'Size;mx'] = Sgam #OK
	1742	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2 #OK
[1459]	1743
	1744	#microstrain derivatives
	1745	if calcControls[phfx+'MustrainType'] == 'isotropic':
[1595]	1746	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
	1747	gamDict[phfx+'Mustrain;i'] = 1.e-6refl[4+im]parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx'] #OK
[1484]	1748	sigDict[phfx+'Mustrain;i'] = 2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])*2/(ateln2parmDict[phfx+'Mustrain;i'])
[1459]	1749	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
	1750	H = np.array(refl[:3])
	1751	P = np.array(calcControls[phfx+'MustrainAxis'])
	1752	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1753	Si = parmDict[phfx+'Mustrain;i']
	1754	Sa = parmDict[phfx+'Mustrain;a']
[1595]	1755	gami = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa
[1459]	1756	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
	1757	Mgam = gami/sqtrm
	1758	dsi = -gamiSicosP2/sqtrm3
	1759	dsa = -gamiSasinP2/sqtrm3
[1484]	1760	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
	1761	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
	1762	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1763	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1764	else: #generalized - P.W. Stephens model OK
[1459]	1765	pwrs = calcControls[phfx+'MuPwrs']
[1479]	1766	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
[1595]	1767	const = 1.e-6parmDict[hfx+'difC']refl[4+im]**3
[1484]	1768	Sum = 0
[1479]	1769	for i,strm in enumerate(Strms):
[1484]	1770	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
	1771	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
	1772	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
	1773	Mgam = const*np.sqrt(Sum)
[1479]	1774	for i in range(len(Strms)):
[1484]	1775	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
	1776	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
[1459]	1777	gamDict[phfx+'Mustrain;mx'] = Mgam
	1778	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
	1779
[1299]	1780	return sigDict,gamDict
	1781
[1597]	1782	def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
[1299]	1783	'Needs a doc string'
[1597]	1784	if im:
	1785	h,k,l,m = refl[:4]
	1786	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
	1787	d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
	1788	else:
	1789	h,k,l = refl[:3]
	1790	d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
[1595]	1791	refl[4+im] = d
[1459]	1792	if 'C' in calcControls[hfx+'histType']:
	1793	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
	1794	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
	1795	if 'Bragg' in calcControls[hfx+'instType']:
	1796	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
	1797	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
	1798	else: #Debye-Scherrer - simple but maybe not right
	1799	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
	1800	elif 'T' in calcControls[hfx+'histType']:
[1475]	1801	pos = parmDict[hfx+'difC']d+parmDict[hfx+'difA']d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
[1459]	1802	#do I need sample position effects - maybe?
[1299]	1803	return pos
	1804
[1597]	1805	def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
[1299]	1806	'Needs a doc string'
	1807	dpr = 180./np.pi
[1597]	1808	if im:
	1809	h,k,l,m = refl[:4]
	1810	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
	1811	dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
	1812	else:
	1813	m = 0
	1814	h,k,l = refl[:3]
	1815	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
[1299]	1816	dst = np.sqrt(dstsq)
[1459]	1817	dsp = 1./dst
	1818	if 'C' in calcControls[hfx+'histType']:
[1595]	1819	pos = refl[5+im]-parmDict[hfx+'Zero']
[1459]	1820	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
	1821	dpdw = const*dst
[1561]	1822	dpdA = np.array([h2,k2,l*2,hk,hl,kl])constwave/(2.0*dst)
[1459]	1823	dpdZ = 1.0
[1597]	1824	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
	1825	2lA[2]+hA[4]+kA[5]])mconstwave/(2.0dst)
	1826	shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
[1459]	1827	if 'Bragg' in calcControls[hfx+'instType']:
[1597]	1828	dpdSh = -4.shftcosd(pos/2.0)
	1829	dpdTr = -shftsind(pos)100.0
	1830	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
[1459]	1831	else: #Debye-Scherrer - simple but maybe not right
[1597]	1832	dpdXd = -shft*cosd(pos)
	1833	dpdYd = -shft*sind(pos)
	1834	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
[1459]	1835	elif 'T' in calcControls[hfx+'histType']:
[1482]	1836	dpdA = -np.array([h2,k2,l*2,hk,hl,kl])parmDict[hfx+'difC']dsp**3/2.
[1459]	1837	dpdZ = 1.0
	1838	dpdDC = dsp
	1839	dpdDA = dsp**2
[1475]	1840	dpdDB = 1./dsp
[1597]	1841	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
	1842	2lA[2]+hA[4]+kA[5]])mparmDict[hfx+'difC']dsp**3/2.
	1843	return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
[1299]	1844
[1595]	1845	def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
[1299]	1846	'Needs a doc string'
	1847	laue = SGData['SGLaue']
	1848	uniq = SGData['SGUniq']
	1849	h,k,l = refl[:3]
	1850	if laue in ['m3','m3m']:
	1851	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
[1595]	1852	refl[4+im]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
[1299]	1853	elif laue in ['6/m','6/mmm','3m1','31m','3']:
	1854	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
	1855	elif laue in ['3R','3mR']:
	1856	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
	1857	elif laue in ['4/m','4/mmm']:
	1858	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
	1859	elif laue in ['mmm']:
	1860	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
	1861	elif laue in ['2/m']:
	1862	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
	1863	if uniq == 'a':
	1864	Dij += parmDict[phfx+'D23']kl
	1865	elif uniq == 'b':
	1866	Dij += parmDict[phfx+'D13']hl
	1867	elif uniq == 'c':
	1868	Dij += parmDict[phfx+'D12']hk
	1869	else:
	1870	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
	1871	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
[1478]	1872	if 'C' in calcControls[hfx+'histType']:
[1595]	1873	return -180.Dijrefl[4+im]*2tand(refl[5+im]/2.0)/np.pi
[1478]	1874	else:
[1595]	1875	return -DijparmDict[hfx+'difC']refl[4+im]**2/2.
[1299]	1876
[1595]	1877	def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
[1299]	1878	'Needs a doc string'
	1879	laue = SGData['SGLaue']
	1880	uniq = SGData['SGUniq']
	1881	h,k,l = refl[:3]
	1882	if laue in ['m3','m3m']:
	1883	dDijDict = {phfx+'D11':h2+k2+l**2,
[1595]	1884	phfx+'eA':refl[4+im]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
[1299]	1885	elif laue in ['6/m','6/mmm','3m1','31m','3']:
	1886	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
	1887	elif laue in ['3R','3mR']:
	1888	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
	1889	elif laue in ['4/m','4/mmm']:
	1890	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
	1891	elif laue in ['mmm']:
	1892	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
	1893	elif laue in ['2/m']:
	1894	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
	1895	if uniq == 'a':
	1896	dDijDict[phfx+'D23'] = k*l
	1897	elif uniq == 'b':
	1898	dDijDict[phfx+'D13'] = h*l
	1899	elif uniq == 'c':
	1900	dDijDict[phfx+'D12'] = h*k
	1901	else:
	1902	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
	1903	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
[1478]	1904	if 'C' in calcControls[hfx+'histType']:
	1905	for item in dDijDict:
[1595]	1906	dDijDict[item] = 180.0refl[4+im]*2tand(refl[5+im]/2.0)/np.pi
[1478]	1907	else:
	1908	for item in dDijDict:
[1595]	1909	dDijDict[item] = -parmDict[hfx+'difC']refl[4+im]**3/2.
[1299]	1910	return dDijDict
	1911
[1559]	1912	def GetDij(phfx,SGData,parmDict):
	1913	HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
	1914	return G2spc.HStrainVals(HSvals,SGData)
	1915
[1299]	1916	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
	1917	'Needs a doc string'
	1918	histoList = Histograms.keys()
	1919	histoList.sort()
	1920	for histogram in histoList:
	1921	if 'PWDR' in histogram[:4]:
	1922	Histogram = Histograms[histogram]
	1923	hId = Histogram['hId']
	1924	hfx = ':%d:'%(hId)
	1925	Limits = calcControls[hfx+'Limits']
[1459]	1926	if 'C' in calcControls[hfx+'histType']:
	1927	shl = max(parmDict[hfx+'SH/L'],0.0005)
	1928	Ka2 = False
	1929	kRatio = 0.0
	1930	if hfx+'Lam1' in parmDict.keys():
	1931	Ka2 = True
	1932	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
	1933	kRatio = parmDict[hfx+'I(L2)/I(L1)']
[1299]	1934	x,y,w,yc,yb,yd = Histogram['Data']
	1935	xB = np.searchsorted(x,Limits[0])
	1936	xF = np.searchsorted(x,Limits[1])
	1937	ymb = np.array(y-yb)
	1938	ymb = np.where(ymb,ymb,1.0)
	1939	ycmb = np.array(yc-yb)
	1940	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
	1941	refLists = Histogram['Reflection Lists']
	1942	for phase in refLists:
[1606]	1943	if phase not in Phases: #skips deleted or renamed phases silently!
	1944	continue
[1299]	1945	Phase = Phases[phase]
[1595]	1946	im = 0
	1947	if Phase['General']['Type'] in ['modulated','magnetic']:
	1948	im = 1
[1299]	1949	pId = Phase['pId']
	1950	phfx = '%d:%d:'%(pId,hId)
	1951	refDict = refLists[phase]
	1952	sumFo = 0.0
	1953	sumdF = 0.0
	1954	sumFosq = 0.0
	1955	sumdFsq = 0.0
[1682]	1956	sumInt = 0.0
[1299]	1957	for refl in refDict['RefList']:
	1958	if 'C' in calcControls[hfx+'histType']:
	1959	yp = np.zeros_like(yb)
[1595]	1960	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
	1961	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
	1962	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
[1299]	1963	iFin2 = iFin
[1498]	1964	if not iBeg+iFin: #peak below low limit - skip peak
	1965	continue
	1966	elif not iBeg-iFin: #peak above high limit - done
	1967	break
	1968	elif iBeg < iFin:
[1595]	1969	yp[iBeg:iFin] = refl[11+im]refl[9+im]G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
[1682]	1970	sumInt += refl[11+im]*refl[9+im]
[1498]	1971	if Ka2:
[1595]	1972	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
	1973	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
[1498]	1974	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
	1975	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
[1798]	1976	if iFin2 > iBeg2:
	1977	yp[iBeg2:iFin2] += refl[11+im]refl[9+im]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2])) #and here
	1978	sumInt += refl[11+im]refl[9+im]kRatio
[1595]	1979	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11+im](1.+kRatio)),0.0))
[1682]	1980
[1299]	1981	elif 'T' in calcControls[hfx+'histType']:
[1459]	1982	yp = np.zeros_like(yb)
[1595]	1983	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
	1984	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
	1985	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
[1493]	1986	if iBeg < iFin:
[1595]	1987	yp[iBeg:iFin] = refl[11+im]refl[9+im]G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin])) #>90% of time spent here
	1988	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
[1682]	1989	sumInt += refl[11+im]*refl[9+im]
[1595]	1990	Fo = np.sqrt(np.abs(refl[8+im]))
	1991	Fc = np.sqrt(np.abs(refl[9]+im))
[1299]	1992	sumFo += Fo
[1595]	1993	sumFosq += refl[8+im]**2
[1299]	1994	sumdF += np.abs(Fo-Fc)
[1595]	1995	sumdFsq += (refl[8+im]-refl[9+im])**2
[1606]	1996	if sumFo:
	1997	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
	1998	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
	1999	else:
	2000	Histogram['Residuals'][phfx+'Rf'] = 100.
	2001	Histogram['Residuals'][phfx+'Rf^2'] = 100.
[1682]	2002	Histogram['Residuals'][phfx+'sumInt'] = sumInt
[1299]	2003	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
	2004	Histogram['Residuals']['hId'] = hId
	2005	elif 'HKLF' in histogram[:4]:
	2006	Histogram = Histograms[histogram]
	2007	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
	2008
	2009	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
	2010	'Needs a doc string'
	2011
[1595]	2012	def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
[1299]	2013	U = parmDict[hfx+'U']
	2014	V = parmDict[hfx+'V']
	2015	W = parmDict[hfx+'W']
	2016	X = parmDict[hfx+'X']
	2017	Y = parmDict[hfx+'Y']
[1595]	2018	tanPos = tand(refl[5+im]/2.0)
	2019	Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
[1299]	2020	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
	2021	sig = max(0.001,sig)
[1595]	2022	gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam #save peak gamma
[1299]	2023	gam = max(0.001,gam)
	2024	return sig,gam
	2025
[1595]	2026	def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
	2027	sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']refl[4+im]*2+ \
	2028	parmDict[hfx+'sig-2']refl[4+im]4+parmDict[hfx+'sig-q']/refl[4+im]*2
	2029	gam = parmDict[hfx+'X']refl[4+im]+parmDict[hfx+'Y']refl[4+im]**2
	2030	Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
[1478]	2031	sig += Ssig
	2032	gam += Sgam
[1459]	2033	return sig,gam
	2034
[1595]	2035	def GetReflAlpBet(refl,im,hfx,parmDict):
	2036	alp = parmDict[hfx+'alpha']/refl[4+im]
	2037	bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]4+parmDict[hfx+'beta-q']/refl[4+im]2
[1459]	2038	return alp,bet
[1559]	2039
[1299]	2040	hId = Histogram['hId']
	2041	hfx = ':%d:'%(hId)
	2042	bakType = calcControls[hfx+'bakType']
[1682]	2043	yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
[1299]	2044	yc = np.zeros_like(yb)
[1461]	2045	cw = np.diff(x)
	2046	cw = np.append(cw,cw[-1])
[1299]	2047
	2048	if 'C' in calcControls[hfx+'histType']:
	2049	shl = max(parmDict[hfx+'SH/L'],0.002)
	2050	Ka2 = False
	2051	if hfx+'Lam1' in parmDict.keys():
	2052	wave = parmDict[hfx+'Lam1']
	2053	Ka2 = True
	2054	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
	2055	kRatio = parmDict[hfx+'I(L2)/I(L1)']
	2056	else:
	2057	wave = parmDict[hfx+'Lam']
	2058	for phase in Histogram['Reflection Lists']:
	2059	refDict = Histogram['Reflection Lists'][phase]
[1606]	2060	if phase not in Phases: #skips deleted or renamed phases silently!
	2061	continue
[1299]	2062	Phase = Phases[phase]
	2063	pId = Phase['pId']
	2064	pfx = '%d::'%(pId)
	2065	phfx = '%d:%d:'%(pId,hId)
	2066	hfx = ':%d:'%(hId)
	2067	SGData = Phase['General']['SGData']
	2068	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
[1595]	2069	im = 0
	2070	if Phase['General']['Type'] in ['modulated','magnetic']:
	2071	SSGData = Phase['General']['SSGData']
	2072	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
	2073	im = 1 #offset in SS reflection list
	2074	#??
[1560]	2075	Dij = GetDij(phfx,SGData,parmDict)
	2076	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
[1299]	2077	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
[1815]	2078	if np.any(np.diag(G)<0.):
	2079	raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
[1299]	2080	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
	2081	Vst = np.sqrt(nl.det(G)) #V*
	2082	if not Phase['General'].get('doPawley'):
	2083	time0 = time.time()
[1618]	2084	if im:
[1630]	2085	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
[1618]	2086	else:
	2087	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
[1299]	2088	# print 'sf calc time: %.3fs'%(time.time()-time0)
	2089	time0 = time.time()
[1495]	2090	badPeak = False
[1299]	2091	for iref,refl in enumerate(refDict['RefList']):
	2092	if 'C' in calcControls[hfx+'histType']:
[1595]	2093	if im:
	2094	h,k,l,m = refl[:4]
	2095	else:
	2096	h,k,l = refl[:3]
[1299]	2097	Uniq = np.inner(refl[:3],SGMT)
[1597]	2098	refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
[1595]	2099	Lorenz = 1./(2.sind(refl[5+im]/2.)2cosd(refl[5+im]/2.)) #Lorentz correction
	2100	refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam
	2101	refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
	2102	refl[11+im] = VstLorenz
[1459]	2103
[1299]	2104	if Phase['General'].get('doPawley'):
	2105	try:
[1595]	2106	if im:
	2107	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
	2108	else:
	2109	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
	2110	refl[9+im] = parmDict[pInd]
[1299]	2111	except KeyError:
	2112	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
	2113	continue
[1595]	2114	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
	2115	iBeg = np.searchsorted(x,refl[5+im]-fmin)
	2116	iFin = np.searchsorted(x,refl[5+im]+fmax)
[1299]	2117	if not iBeg+iFin: #peak below low limit - skip peak
	2118	continue
	2119	elif not iBeg-iFin: #peak above high limit - done
	2120	break
[1494]	2121	elif iBeg > iFin: #bad peak coeff - skip
[1495]	2122	badPeak = True
[1494]	2123	continue
[1595]	2124	yc[iBeg:iFin] += refl[11+im]refl[9+im]G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
[1299]	2125	if Ka2:
[1595]	2126	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
	2127	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
[1299]	2128	iBeg = np.searchsorted(x,pos2-fmin)
	2129	iFin = np.searchsorted(x,pos2+fmax)
	2130	if not iBeg+iFin: #peak below low limit - skip peak
	2131	continue
	2132	elif not iBeg-iFin: #peak above high limit - done
	2133	return yc,yb
[1494]	2134	elif iBeg > iFin: #bad peak coeff - skip
	2135	continue
[1595]	2136	yc[iBeg:iFin] += refl[11+im]refl[9+im]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #and here
[1299]	2137	elif 'T' in calcControls[hfx+'histType']:
[1459]	2138	h,k,l = refl[:3]
	2139	Uniq = np.inner(refl[:3],SGMT)
[1597]	2140	refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
[1766]	2141	Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)refl[4+im]*4 #TOF Lorentz correction
[1595]	2142	# refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift
	2143	refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict) #peak sig & gam
	2144	refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
	2145	refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
	2146	refl[11+im] = VstLorenz
[1459]	2147	if Phase['General'].get('doPawley'):
	2148	try:
[1595]	2149	if im:
	2150	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
	2151	else:
	2152	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
	2153	refl[9+im] = parmDict[pInd]
[1459]	2154	except KeyError:
	2155	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
	2156	continue
[1595]	2157	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
	2158	iBeg = np.searchsorted(x,refl[5+im]-fmin)
	2159	iFin = np.searchsorted(x,refl[5+im]+fmax)
[1459]	2160	if not iBeg+iFin: #peak below low limit - skip peak
	2161	continue
	2162	elif not iBeg-iFin: #peak above high limit - done
	2163	break
[1494]	2164	elif iBeg > iFin: #bad peak coeff - skip
[1495]	2165	badPeak = True
[1494]	2166	continue
[1595]	2167	yc[iBeg:iFin] += refl[11+im]refl[9+im]G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
[1299]	2168	# print 'profile calc time: %.3fs'%(time.time()-time0)
[1495]	2169	if badPeak:
	2170	print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
[1299]	2171	return yc,yb
	2172
	2173	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
	2174	'Needs a doc string'
	2175
	2176	def cellVaryDerv(pfx,SGData,dpdA):
	2177	if SGData['SGLaue'] in ['-1',]:
	2178	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
	2179	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
	2180	elif SGData['SGLaue'] in ['2/m',]:
	2181	if SGData['SGUniq'] == 'a':
[1762]	2182	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
[1299]	2183	elif SGData['SGUniq'] == 'b':
	2184	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
	2185	else:
[1762]	2186	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
[1299]	2187	elif SGData['SGLaue'] in ['mmm',]:
	2188	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
	2189	elif SGData['SGLaue'] in ['4/m','4/mmm']:
	2190	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
	2191	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
	2192	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
	2193	elif SGData['SGLaue'] in ['3R', '3mR']:
	2194	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
	2195	elif SGData['SGLaue'] in ['m3m','m3']:
	2196	return [[pfx+'A0',dpdA[0]]]
	2197
	2198	# create a list of dependent variables and set up a dictionary to hold their derivatives
	2199	dependentVars = G2mv.GetDependentVars()
	2200	depDerivDict = {}
	2201	for j in dependentVars:
	2202	depDerivDict[j] = np.zeros(shape=(len(x)))
	2203	#print 'dependent vars',dependentVars
	2204	lenX = len(x)
	2205	hId = Histogram['hId']
	2206	hfx = ':%d:'%(hId)
	2207	bakType = calcControls[hfx+'bakType']
	2208	dMdv = np.zeros(shape=(len(varylist),len(x)))
[1474]	2209	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
[1496]	2210	if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
	2211	bBpos =varylist.index(hfx+'Back;0')
[1299]	2212	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
	2213	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
	2214	for name in varylist:
	2215	if 'Debye' in name:
[1496]	2216	id = int(name.split(';')[-1])
	2217	parm = name[:int(name.rindex(';'))]
[1299]	2218	ip = names.index(parm)
	2219	dMdv[varylist.index(name)] = dMddb[3*id+ip]
	2220	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
	2221	for name in varylist:
	2222	if 'BkPk' in name:
	2223	parm,id = name.split(';')
	2224	id = int(id)
	2225	if parm in names:
	2226	ip = names.index(parm)
	2227	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
	2228	cw = np.diff(x)
	2229	cw = np.append(cw,cw[-1])
[1459]	2230	Ka2 = False #also for TOF!
[1299]	2231	if 'C' in calcControls[hfx+'histType']:
	2232	shl = max(parmDict[hfx+'SH/L'],0.002)
	2233	if hfx+'Lam1' in parmDict.keys():
	2234	wave = parmDict[hfx+'Lam1']
	2235	Ka2 = True
	2236	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
	2237	kRatio = parmDict[hfx+'I(L2)/I(L1)']
	2238	else:
	2239	wave = parmDict[hfx+'Lam']
	2240	for phase in Histogram['Reflection Lists']:
	2241	refDict = Histogram['Reflection Lists'][phase]
[1606]	2242	if phase not in Phases: #skips deleted or renamed phases silently!
	2243	continue
[1299]	2244	Phase = Phases[phase]
	2245	SGData = Phase['General']['SGData']
	2246	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
[1595]	2247	im = 0
	2248	if Phase['General']['Type'] in ['modulated','magnetic']:
	2249	SSGData = Phase['General']['SSGData']
	2250	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
	2251	im = 1 #offset in SS reflection list
	2252	#??
[1299]	2253	pId = Phase['pId']
	2254	pfx = '%d::'%(pId)
	2255	phfx = '%d:%d:'%(pId,hId)
[1560]	2256	Dij = GetDij(phfx,SGData,parmDict)
	2257	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
[1299]	2258	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	2259	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
	2260	if not Phase['General'].get('doPawley'):
	2261	time0 = time.time()
[1618]	2262	if im:
[1630]	2263	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
[1618]	2264	else:
	2265	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
[1299]	2266	# print 'sf-derv time %.3fs'%(time.time()-time0)
	2267	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
	2268	time0 = time.time()
	2269	for iref,refl in enumerate(refDict['RefList']):
[1597]	2270	if im:
	2271	h,k,l,m = refl[:4]
	2272	else:
	2273	h,k,l = refl[:3]
[1459]	2274	Uniq = np.inner(refl[:3],SGMT)
[1474]	2275	if 'T' in calcControls[hfx+'histType']:
[1595]	2276	wave = refl[14+im]
	2277	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
[1299]	2278	if 'C' in calcControls[hfx+'histType']: #CW powder
[1595]	2279	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
[1459]	2280	else: #'T'OF
[1595]	2281	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
	2282	iBeg = np.searchsorted(x,refl[5+im]-fmin)
	2283	iFin = np.searchsorted(x,refl[5+im]+fmax)
[1459]	2284	if not iBeg+iFin: #peak below low limit - skip peak
	2285	continue
	2286	elif not iBeg-iFin: #peak above high limit - done
	2287	break
[1595]	2288	pos = refl[5+im]
[1459]	2289	if 'C' in calcControls[hfx+'histType']:
[1299]	2290	tanth = tand(pos/2.0)
	2291	costh = cosd(pos/2.0)
	2292	lenBF = iFin-iBeg
	2293	dMdpk = np.zeros(shape=(6,lenBF))
[1595]	2294	dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
[1299]	2295	for i in range(5):
[1595]	2296	dMdpk[i] += 100.cw[iBeg:iFin]refl[11+im]refl[9+im]dMdipk[i]
[1299]	2297	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
	2298	if Ka2:
[1595]	2299	pos2 = refl[5+im]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
[1299]	2300	iBeg2 = np.searchsorted(x,pos2-fmin)
	2301	iFin2 = np.searchsorted(x,pos2+fmax)
	2302	if iBeg2-iFin2:
	2303	lenBF2 = iFin2-iBeg2
	2304	dMdpk2 = np.zeros(shape=(6,lenBF2))
[1595]	2305	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
[1299]	2306	for i in range(5):
[1595]	2307	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11+im]refl[9+im]kRatio*dMdipk2[i]
	2308	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11+im]*dMdipk2[0]
[1299]	2309	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
[1459]	2310	else: #'T'OF
	2311	lenBF = iFin-iBeg
[1494]	2312	if lenBF < 0: #bad peak coeff
	2313	break
[1459]	2314	dMdpk = np.zeros(shape=(6,lenBF))
[1595]	2315	dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
[1460]	2316	for i in range(6):
[1595]	2317	dMdpk[i] += refl[11+im]refl[9+im]dMdipk[i] #cw[iBeg:iFin]*
[1459]	2318	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
	2319	if Phase['General'].get('doPawley'):
	2320	dMdpw = np.zeros(len(x))
	2321	try:
[1595]	2322	if im:
	2323	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
	2324	else:
	2325	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
[1459]	2326	idx = varylist.index(pIdx)
[1595]	2327	dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
[1459]	2328	if Ka2: #not for TOF either
[1595]	2329	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
[1459]	2330	dMdv[idx] = dMdpw
	2331	except: # ValueError:
	2332	pass
	2333	if 'C' in calcControls[hfx+'histType']:
[1597]	2334	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
[1299]	2335	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
	2336	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
	2337	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
	2338	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
	2339	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
	2340	hfx+'DisplaceY':[dpdY,'pos'],}
	2341	if 'Bragg' in calcControls[hfx+'instType']:
	2342	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
	2343	hfx+'SurfRoughB':[dFdAb[1],'int'],})
	2344	else:
	2345	names.update({hfx+'Absorption':[dFdAb,'int'],})
[1459]	2346	else: #'T'OF
[1597]	2347	dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
[1459]	2348	names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
	2349	hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
[1595]	2350	hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
	2351	hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
	2352	hfx+'beta-q':[1./refl[4+im]2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]2,'sig'],
	2353	hfx+'sig-2':[refl[4+im]4,'sig'],hfx+'sig-q':[1./refl[4+im]2,'sig'],
[1465]	2354	hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
[1459]	2355	for name in names:
	2356	item = names[name]
	2357	if name in varylist:
	2358	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
[1798]	2359	if Ka2 and iFin2-iBeg2:
[1459]	2360	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
	2361	elif name in dependentVars:
	2362	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
[1798]	2363	if Ka2 and iFin2-iBeg2:
[1459]	2364	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
	2365	for iPO in dIdPO:
	2366	if iPO in varylist:
	2367	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
[1798]	2368	if Ka2 and iFin2-iBeg2:
[1459]	2369	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
	2370	elif iPO in dependentVars:
	2371	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
[1798]	2372	if Ka2 and iFin2-iBeg2:
[1459]	2373	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
	2374	for i,name in enumerate(['omega','chi','phi']):
	2375	aname = pfx+'SH '+name
	2376	if aname in varylist:
	2377	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
[1798]	2378	if Ka2 and iFin2-iBeg2:
[1459]	2379	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
	2380	elif aname in dependentVars:
	2381	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
[1798]	2382	if Ka2 and iFin2-iBeg2:
[1459]	2383	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
	2384	for iSH in dFdODF:
	2385	if iSH in varylist:
	2386	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
[1798]	2387	if Ka2 and iFin2-iBeg2:
[1459]	2388	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
	2389	elif iSH in dependentVars:
	2390	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
[1798]	2391	if Ka2 and iFin2-iBeg2:
[1459]	2392	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
	2393	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
	2394	for name,dpdA in cellDervNames:
	2395	if name in varylist:
	2396	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
[1798]	2397	if Ka2 and iFin2-iBeg2:
[1459]	2398	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
	2399	elif name in dependentVars:
	2400	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
[1798]	2401	if Ka2 and iFin2-iBeg2:
[1459]	2402	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
[1595]	2403	dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
[1459]	2404	for name in dDijDict:
	2405	if name in varylist:
	2406	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
[1798]	2407	if Ka2 and iFin2-iBeg2:
[1459]	2408	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
	2409	elif name in dependentVars:
	2410	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
[1798]	2411	if Ka2 and iFin2-iBeg2:
[1459]	2412	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
[1597]	2413	for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
	2414	if name in varylist:
	2415	dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
[1798]	2416	if Ka2 and iFin2-iBeg2:
[1597]	2417	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
	2418	elif name in dependentVars:
	2419	depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
[1798]	2420	if Ka2 and iFin2-iBeg2:
[1597]	2421	depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
[1459]	2422	if 'C' in calcControls[hfx+'histType']:
[1595]	2423	sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
[1459]	2424	else: #'T'OF
[1595]	2425	sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
[1459]	2426	for name in gamDict:
	2427	if name in varylist:
	2428	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
[1798]	2429	if Ka2 and iFin2-iBeg2:
[1459]	2430	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
	2431	elif name in dependentVars:
	2432	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
[1798]	2433	if Ka2 and iFin2-iBeg2:
[1459]	2434	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
	2435	for name in sigDict:
	2436	if name in varylist:
	2437	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
[1798]	2438	if Ka2 and iFin2-iBeg2:
[1459]	2439	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
	2440	elif name in dependentVars:
	2441	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
[1798]	2442	if Ka2 and iFin2-iBeg2:
[1459]	2443	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
	2444	for name in ['BabA','BabU']:
[1595]	2445	if refl[9+im]:
[1459]	2446	if phfx+name in varylist:
[1595]	2447	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
[1798]	2448	if Ka2 and iFin2-iBeg2:
[1595]	2449	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
[1459]	2450	elif phfx+name in dependentVars:
[1595]	2451	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
[1798]	2452	if Ka2 and iFin2-iBeg2:
[1595]	2453	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
[1299]	2454	if not Phase['General'].get('doPawley'):
	2455	#do atom derivatives - for RB,F,X & U so far
[1595]	2456	corr = dervDict['int']/refl[9+im]
[1798]	2457	if Ka2 and iFin2-iBeg2:
[1595]	2458	corr2 = dervDict2['int']/refl[9+im]
[1299]	2459	for name in varylist+dependentVars:
	2460	if '::RBV;' in name:
	2461	pass
	2462	else:
	2463	try:
	2464	aname = name.split(pfx)[1][:2]
	2465	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
	2466	except IndexError:
	2467	continue
	2468	if name in varylist:
	2469	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
[1798]	2470	if Ka2 and iFin2-iBeg2:
[1299]	2471	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
	2472	elif name in dependentVars:
	2473	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
[1798]	2474	if Ka2 and iFin2-iBeg2:
[1299]	2475	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
[1459]	2476	# print 'profile derv time: %.3fs'%(time.time()-time0)
[1299]	2477	# now process derivatives in constraints
	2478	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
	2479	return dMdv
[1456]	2480
[1787]	2481	def UserRejectHKL(ref,im,userReject):
	2482	if ref[5+im]/ref[6+im] < userReject['minF/sig']:
	2483	return False
	2484	elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
	2485	return False
	2486	elif ref[11+im] < userReject['MinExt']:
	2487	return False
	2488	elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
	2489	return False
	2490	return True
	2491
[1456]	2492	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
[1459]	2493	'''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
[1456]	2494	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
	2495	are returned to either dervRefine or HessRefine.
[1299]	2496
[1456]	2497	:returns:
	2498	'''
	2499	nobs = Histogram['Residuals']['Nobs']
	2500	hId = Histogram['hId']
	2501	hfx = ':%d:'%(hId)
	2502	pfx = '%d::'%(Phase['pId'])
	2503	phfx = '%d:%d:'%(Phase['pId'],hId)
	2504	SGData = Phase['General']['SGData']
[1595]	2505	im = 0
	2506	if Phase['General']['Type'] in ['modulated','magnetic']:
	2507	SSGData = Phase['General']['SSGData']
	2508	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
	2509	im = 1 #offset in SS reflection list
[1456]	2510	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	2511	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	2512	refDict = Histogram['Data']
[1618]	2513	if im:
[1630]	2514	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
[1618]	2515	else:
	2516	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
[1456]	2517	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
	2518	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
	2519	dependentVars = G2mv.GetDependentVars()
	2520	depDerivDict = {}
	2521	for j in dependentVars:
	2522	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
	2523	wdf = np.zeros(len(refDict['RefList']))
	2524	if calcControls['F**2']:
	2525	for iref,ref in enumerate(refDict['RefList']):
[1595]	2526	if ref[6+im] > 0:
[1857]	2527	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
[1595]	2528	w = 1.0/ref[6+im]
[1787]	2529	if ref[3+im] > 0:
[1595]	2530	wdf[iref] = w*(ref[5+im]-ref[7+im])
[1456]	2531	for j,var in enumerate(varylist):
	2532	if var in dFdvDict:
[1595]	2533	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2534	for var in dependentVars:
	2535	if var in dFdvDict:
[1595]	2536	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2537	if phfx+'Scale' in varylist:
[1857]	2538	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9+im]ref[11+im] #OK
[1456]	2539	elif phfx+'Scale' in dependentVars:
[1857]	2540	depDerivDict[phfx+'Scale'][iref] = wref[9+im]ref[11+im] #OK
[1456]	2541	for item in ['Ep','Es','Eg']:
[1493]	2542	if phfx+item in varylist and phfx+item in dervDict:
[1595]	2543	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
[1493]	2544	elif phfx+item in dependentVars and phfx+item in dervDict:
[1595]	2545	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
[1456]	2546	for item in ['BabA','BabU']:
	2547	if phfx+item in varylist:
[1595]	2548	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2549	elif phfx+item in dependentVars:
[1595]	2550	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
[1857]	2551	else: #F refinement
[1456]	2552	for iref,ref in enumerate(refDict['RefList']):
[1595]	2553	if ref[5+im] > 0.:
[1857]	2554	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
[1595]	2555	Fo = np.sqrt(ref[5+im])
	2556	Fc = np.sqrt(ref[7+im])
	2557	w = 1.0/ref[6+im]
[1787]	2558	if ref[3+im] > 0:
[1864]	2559	wdf[iref] = 2.0Fcw*(Fo-Fc)
[1456]	2560	for j,var in enumerate(varylist):
	2561	if var in dFdvDict:
[1595]	2562	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2563	for var in dependentVars:
	2564	if var in dFdvDict:
[1595]	2565	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2566	if phfx+'Scale' in varylist:
[1864]	2567	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9+im]ref[11+im] #OK
[1456]	2568	elif phfx+'Scale' in dependentVars:
[1864]	2569	depDerivDict[phfx+'Scale'][iref] = wref[9+im]ref[11+im] #OK
	2570	for item in ['Ep','Es','Eg']: #OK!
[1493]	2571	if phfx+item in varylist and phfx+item in dervDict:
[1862]	2572	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]
[1493]	2573	elif phfx+item in dependentVars and phfx+item in dervDict:
[1862]	2574	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
[1456]	2575	for item in ['BabA','BabU']:
	2576	if phfx+item in varylist:
[1595]	2577	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2578	elif phfx+item in dependentVars:
[1595]	2579	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
[1456]	2580	return dMdvh,depDerivDict,wdf
	2581
	2582
[1299]	2583	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
	2584	'''Loop over histograms and compute derivatives of the fitting
	2585	model (M) with respect to all parameters. Results are returned in
	2586	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
	2587	n is the number of parameters and m is the number of data
	2588	points. This can exceed memory when m gets large. This routine is
	2589	used when refinement derivatives are selected as "analtytic
	2590	Jacobian" in Controls.
	2591
	2592	:returns: Jacobian numpy.array dMdv for all histograms concatinated
	2593	'''
	2594	parmDict.update(zip(varylist,values))
	2595	G2mv.Dict2Map(parmDict,varylist)
	2596	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	2597	nvar = len(varylist)
	2598	dMdv = np.empty(0)
	2599	histoList = Histograms.keys()
	2600	histoList.sort()
	2601	for histogram in histoList:
	2602	if 'PWDR' in histogram[:4]:
	2603	Histogram = Histograms[histogram]
	2604	hId = Histogram['hId']
	2605	hfx = ':%d:'%(hId)
	2606	wtFactor = calcControls[hfx+'wtFactor']
	2607	Limits = calcControls[hfx+'Limits']
	2608	x,y,w,yc,yb,yd = Histogram['Data']
	2609	xB = np.searchsorted(x,Limits[0])
	2610	xF = np.searchsorted(x,Limits[1])
[1489]	2611	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
[1299]	2612	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
	2613	elif 'HKLF' in histogram[:4]:
	2614	Histogram = Histograms[histogram]
	2615	phase = Histogram['Reflection Lists']
	2616	Phase = Phases[phase]
[1456]	2617	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
	2618	hfx = ':%d:'%(Histogram['hId'])
[1299]	2619	wtFactor = calcControls[hfx+'wtFactor']
	2620	# now process derivatives in constraints
	2621	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
	2622	else:
	2623	continue #skip non-histogram entries
	2624	if len(dMdv):
	2625	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
	2626	else:
	2627	dMdv = np.sqrt(wtFactor)*dMdvh
	2628
[1775]	2629	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
[1299]	2630	if np.any(pVals):
[1775]	2631	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
[1299]	2632	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
	2633
	2634	return dMdv
	2635
	2636	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
	2637	'''Loop over histograms and compute derivatives of the fitting
	2638	model (M) with respect to all parameters. For each histogram, the
	2639	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
	2640	number of parameters and m is the number of data points *in the
	2641	histogram*. The (n by n) Hessian is computed from each Jacobian
	2642	and it is returned. This routine is used when refinement
	2643	derivatives are selected as "analtytic Hessian" in Controls.
	2644
	2645	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
	2646	'''
	2647	parmDict.update(zip(varylist,values))
	2648	G2mv.Dict2Map(parmDict,varylist)
	2649	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	2650	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
	2651	nvar = len(varylist)
	2652	Hess = np.empty(0)
	2653	histoList = Histograms.keys()
	2654	histoList.sort()
	2655	for histogram in histoList:
	2656	if 'PWDR' in histogram[:4]:
	2657	Histogram = Histograms[histogram]
	2658	hId = Histogram['hId']
	2659	hfx = ':%d:'%(hId)
	2660	wtFactor = calcControls[hfx+'wtFactor']
	2661	Limits = calcControls[hfx+'Limits']
	2662	x,y,w,yc,yb,yd = Histogram['Data']
	2663	W = wtFactor*w
	2664	dy = y-yc
	2665	xB = np.searchsorted(x,Limits[0])
	2666	xF = np.searchsorted(x,Limits[1])
	2667	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
	2668	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
	2669	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
	2670	Dy = dy[xB:xF][np.newaxis,:]
	2671	dMdvh *= Wt
	2672	if dlg:
[1591]	2673	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
[1299]	2674	if len(Hess):
	2675	Hess += np.inner(dMdvh,dMdvh)
	2676	dMdvh = WtDy
	2677	Vec += np.sum(dMdvh,axis=1)
	2678	else:
	2679	Hess = np.inner(dMdvh,dMdvh)
	2680	dMdvh = WtDy
	2681	Vec = np.sum(dMdvh,axis=1)
	2682	elif 'HKLF' in histogram[:4]:
	2683	Histogram = Histograms[histogram]
	2684	phase = Histogram['Reflection Lists']
	2685	Phase = Phases[phase]
[1456]	2686	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
[1299]	2687	hId = Histogram['hId']
[1456]	2688	hfx = ':%d:'%(Histogram['hId'])
[1299]	2689	wtFactor = calcControls[hfx+'wtFactor']
	2690	# now process derivatives in constraints
	2691	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
[1453]	2692	# print 'matrix build time: %.3f'%(time.time()-time0)
[1299]	2693
	2694	if dlg:
	2695	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
	2696	if len(Hess):
	2697	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
	2698	Hess += wtFactor*np.inner(dMdvh,dMdvh)
	2699	else:
	2700	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
	2701	Hess = wtFactor*np.inner(dMdvh,dMdvh)
	2702	else:
	2703	continue #skip non-histogram entries
[1775]	2704	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
[1299]	2705	if np.any(pVals):
[1775]	2706	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
[1299]	2707	Vec += np.sum(dpdvpWtpVals,axis=1)
	2708	Hess += np.inner(dpdv*pWt,dpdv)
	2709	return Vec,Hess
	2710
[1676]	2711	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):
	2712	'''Computes the point-by-point discrepancies between every data point in every histogram
	2713	and the observed value
	2714
	2715	:returns: an np array of differences between observed and computed diffraction values.
	2716	'''
[1322]	2717	Values2Dict(parmDict, varylist, values)
[1299]	2718	G2mv.Dict2Map(parmDict,varylist)
	2719	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	2720	M = np.empty(0)
	2721	SumwYo = 0
	2722	Nobs = 0
[1782]	2723	Nrej = 0
[1800]	2724	Next = 0
[1299]	2725	ApplyRBModels(parmDict,Phases,rigidbodyDict)
	2726	histoList = Histograms.keys()
	2727	histoList.sort()
	2728	for histogram in histoList:
	2729	if 'PWDR' in histogram[:4]:
	2730	Histogram = Histograms[histogram]
	2731	hId = Histogram['hId']
	2732	hfx = ':%d:'%(hId)
	2733	wtFactor = calcControls[hfx+'wtFactor']
	2734	Limits = calcControls[hfx+'Limits']
	2735	x,y,w,yc,yb,yd = Histogram['Data']
	2736	yc *= 0.0 #zero full calcd profiles
	2737	yb *= 0.0
	2738	yd *= 0.0
	2739	xB = np.searchsorted(x,Limits[0])
	2740	xF = np.searchsorted(x,Limits[1])
	2741	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
	2742	varylist,Histogram,Phases,calcControls,pawleyLookup)
	2743	yc[xB:xF] += yb[xB:xF]
	2744	if not np.any(y): #fill dummy data
	2745	rv = st.poisson(yc[xB:xF])
	2746	y[xB:xF] = rv.rvs()
	2747	Z = np.ones_like(yc[xB:xF])
	2748	Z[1::2] *= -1
	2749	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
	2750	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
	2751	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
	2752	W = wtFactor*w
	2753	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
	2754	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
	2755	Nobs += Histogram['Residuals']['Nobs']
	2756	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
	2757	SumwYo += Histogram['Residuals']['sumwYo']
	2758	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
	2759	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
	2760	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
	2761	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
	2762	sumYB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yc[xB:xF])*ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF],0.))
	2763	sumwYB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](ma.abs(y[xB:xF]-yc[xB:xF])ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF])**2,0.))
	2764	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
	2765	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
	2766	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
	2767	if dlg:
	2768	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
	2769	M = np.concatenate((M,wdy))
	2770	#end of PWDR processing
	2771	elif 'HKLF' in histogram[:4]:
	2772	Histogram = Histograms[histogram]
	2773	Histogram['Residuals'] = {}
	2774	phase = Histogram['Reflection Lists']
	2775	Phase = Phases[phase]
	2776	hId = Histogram['hId']
	2777	hfx = ':%d:'%(hId)
	2778	wtFactor = calcControls[hfx+'wtFactor']
	2779	pfx = '%d::'%(Phase['pId'])
	2780	phfx = '%d:%d:'%(Phase['pId'],hId)
	2781	SGData = Phase['General']['SGData']
[1595]	2782	im = 0
	2783	if Phase['General']['Type'] in ['modulated','magnetic']:
	2784	SSGData = Phase['General']['SSGData']
	2785	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
	2786	im = 1 #offset in SS reflection list
[1299]	2787	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	2788	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	2789	refDict = Histogram['Data']
	2790	time0 = time.time()
[1613]	2791	if im:
[1630]	2792	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
[1613]	2793	else:
[1618]	2794	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
[1453]	2795	# print 'sf-calc time: %.3f'%(time.time()-time0)
[1299]	2796	df = np.zeros(len(refDict['RefList']))
	2797	sumwYo = 0
	2798	sumFo = 0
	2799	sumFo2 = 0
	2800	sumdF = 0
	2801	sumdF2 = 0
	2802	nobs = 0
[1782]	2803	nrej = 0
[1800]	2804	next = 0
[1299]	2805	if calcControls['F**2']:
	2806	for i,ref in enumerate(refDict['RefList']):
[1595]	2807	if ref[6+im] > 0:
[1604]	2808	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
[1784]	2809	w = 1.0/ref[6+im] # 1/sig(F^2)
[1886]	2810	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
[1595]	2811	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
[1797]	2812	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
[1787]	2813	ref[3+im] = abs(ref[3+im]) #mark as allowed
[1595]	2814	Fo = np.sqrt(ref[5+im])
[1299]	2815	sumFo += Fo
[1595]	2816	sumFo2 += ref[5+im]
	2817	sumdF += abs(Fo-np.sqrt(ref[7+im]))
	2818	sumdF2 += abs(ref[5+im]-ref[7+im])
[1299]	2819	nobs += 1
[1595]	2820	df[i] = -w*(ref[5+im]-ref[7+im])
[1784]	2821	sumwYo += (wref[5+im])2 #wFo^2
[1782]	2822	else:
[1800]	2823	if ref[3+im]:
	2824	ref[3+im] = -abs(ref[3+im]) #mark as rejected
	2825	nrej += 1
	2826	else: #sp.gp.extinct
	2827	next += 1
[1299]	2828	else:
	2829	for i,ref in enumerate(refDict['RefList']):
[1595]	2830	if ref[5+im] > 0.:
	2831	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
[1886]	2832	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
[1595]	2833	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
	2834	Fo = np.sqrt(ref[5+im])
	2835	Fc = np.sqrt(ref[7+im])
[1864]	2836	w = 2.0*Fo/ref[6+im] # 1/sig(F)?
[1797]	2837	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
[1787]	2838	ref[3+im] = abs(ref[3+im]) #mark as allowed
[1299]	2839	sumFo += Fo
[1595]	2840	sumFo2 += ref[5+im]
[1299]	2841	sumdF += abs(Fo-Fc)
[1595]	2842	sumdF2 += abs(ref[5+im]-ref[7+im])
[1299]	2843	nobs += 1
	2844	df[i] = -w*(Fo-Fc)
	2845	sumwYo += (wFo)*2
[1782]	2846	else:
[1800]	2847	if ref[3+im]:
	2848	ref[3+im] = -abs(ref[3+im]) #mark as rejected
	2849	nrej += 1
	2850	else: #sp.gp.extinct
	2851	next += 1
[1299]	2852	Histogram['Residuals']['Nobs'] = nobs
	2853	Histogram['Residuals']['sumwYo'] = sumwYo
	2854	SumwYo += sumwYo
[1784]	2855	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/sumwYo)100.)
[1299]	2856	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
	2857	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
	2858	Histogram['Residuals'][phfx+'Nref'] = nobs
[1782]	2859	Histogram['Residuals'][phfx+'Nrej'] = nrej
[1800]	2860	Histogram['Residuals'][phfx+'Next'] = next
[1299]	2861	Nobs += nobs
[1782]	2862	Nrej += nrej
[1800]	2863	Next += next
[1299]	2864	if dlg:
	2865	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
	2866	M = np.concatenate((M,wtFactor*df))
	2867	# end of HKLF processing
	2868	Histograms['sumwYo'] = SumwYo
	2869	Histograms['Nobs'] = Nobs
[1782]	2870	Histograms['Nrej'] = Nrej
[1800]	2871	Histograms['Next'] = Next
[1299]	2872	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
	2873	if dlg:
	2874	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
	2875	if not GoOn:
	2876	parmDict['saved values'] = values
	2877	dlg.Destroy()
[1813]	2878	raise G2obj.G2Exception('User abort') #Abort!!
[1775]	2879	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
[1365]	2880	if len(pVals):
[1299]	2881	pSum = np.sum(pWtpVals*2)
	2882	for name in pWsum:
[1775]	2883	if pWsum[name]:
[1365]	2884	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
[1359]	2885	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
[1299]	2886	Nobs += len(pVals)
	2887	M = np.concatenate((M,np.sqrt(pWt)*pVals))
	2888	return M
	2889

Note: See TracBrowser for help on using the repository browser.

Download in other formats: