Context Navigation

source: trunk/GSASIIstrMath.py @ 1088

Last change on this file since 1088 was 1088, checked in by vondreele, 12 years ago
small tweak to data simulation - give it a "signature" of fakeness
Property svn:keywords set to `Date Author Revision URL Id`
File size: 101.0 KB

Rev	Line
[926]	1	# -- coding: utf-8 --
[939]	2	'''
	3	GSASIIstrMath - structure math routines
	4	-----------------------------------------
	5	'''
[926]	6	########### SVN repository information ###################
[1077]	7	# $Date: 2013-10-07 18:01:01 +0000 (Mon, 07 Oct 2013) $
	8	# $Author: vondreele $
	9	# $Revision: 1088 $
	10	# $URL: trunk/GSASIIstrMath.py $
	11	# $Id: GSASIIstrMath.py 1088 2013-10-07 18:01:01Z vondreele $
[926]	12	########### SVN repository information ###################
	13	import time
	14	import math
	15	import copy
	16	import numpy as np
	17	import numpy.ma as ma
	18	import numpy.linalg as nl
	19	import scipy.optimize as so
[1087]	20	import scipy.stats as st
[926]	21	import GSASIIpath
[1077]	22	GSASIIpath.SetVersionNumber("$Revision: 1088 $")
[926]	23	import GSASIIElem as G2el
	24	import GSASIIlattice as G2lat
	25	import GSASIIspc as G2spc
	26	import GSASIIpwd as G2pwd
	27	import GSASIImapvars as G2mv
	28	import GSASIImath as G2mth
	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
	37	ateln2 = 8.0*math.log(2.0)
	38
	39	################################################################################
	40	##### Rigid Body Models
	41	################################################################################
	42
	43	def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
	44	''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
	45	current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
	46	'''
	47	atxIds = ['Ax:','Ay:','Az:']
	48	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
	49	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
	50	if not RBIds['Vector'] and not RBIds['Residue']:
	51	return
	52	VRBIds = RBIds['Vector']
	53	RRBIds = RBIds['Residue']
	54	if Update:
	55	RBData = rigidbodyDict
	56	else:
	57	RBData = copy.deepcopy(rigidbodyDict) # don't mess with original!
	58	if RBIds['Vector']: # first update the vector magnitudes
	59	VRBData = RBData['Vector']
	60	for i,rbId in enumerate(VRBIds):
[961]	61	if VRBData[rbId]['useCount']:
	62	for j in range(len(VRBData[rbId]['VectMag'])):
	63	name = '::RBV;'+str(j)+':'+str(i)
	64	VRBData[rbId]['VectMag'][j] = parmDict[name]
[926]	65	for phase in Phases:
	66	Phase = Phases[phase]
	67	General = Phase['General']
	68	cell = General['Cell'][1:7]
	69	Amat,Bmat = G2lat.cell2AB(cell)
	70	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
	71	pfx = str(Phase['pId'])+'::'
	72	if Update:
	73	RBModels = Phase['RBModels']
	74	else:
	75	RBModels = copy.deepcopy(Phase['RBModels']) # again don't mess with original!
	76	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
	77	jrb = VRBIds.index(RBObj['RBId'])
	78	rbsx = str(irb)+':'+str(jrb)
	79	for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
	80	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
	81	for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
	82	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
[953]	83	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
[926]	84	TLS = RBObj['ThermalMotion']
	85	if 'T' in TLS[0]:
	86	for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
	87	TLS[1][i] = parmDict[pfx+pt+rbsx]
	88	if 'L' in TLS[0]:
	89	for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
	90	TLS[1][i+6] = parmDict[pfx+pt+rbsx]
	91	if 'S' in TLS[0]:
	92	for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
	93	TLS[1][i+12] = parmDict[pfx+pt+rbsx]
	94	if 'U' in TLS[0]:
	95	TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
	96	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
	97	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
	98	for i,x in enumerate(XYZ):
	99	atId = RBObj['Ids'][i]
	100	for j in [0,1,2]:
	101	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
	102	if UIJ[i][0] == 'A':
	103	for j in range(6):
	104	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
	105	elif UIJ[i][0] == 'I':
	106	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
	107
	108	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
	109	jrb = RRBIds.index(RBObj['RBId'])
	110	rbsx = str(irb)+':'+str(jrb)
	111	for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
	112	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
	113	for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
	114	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
[953]	115	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
[926]	116	TLS = RBObj['ThermalMotion']
	117	if 'T' in TLS[0]:
	118	for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
	119	RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
	120	if 'L' in TLS[0]:
	121	for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
	122	RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
	123	if 'S' in TLS[0]:
	124	for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
	125	RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
	126	if 'U' in TLS[0]:
	127	RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
	128	for itors,tors in enumerate(RBObj['Torsions']):
	129	tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
	130	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
	131	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
	132	for i,x in enumerate(XYZ):
	133	atId = RBObj['Ids'][i]
	134	for j in [0,1,2]:
	135	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
	136	if UIJ[i][0] == 'A':
	137	for j in range(6):
	138	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
	139	elif UIJ[i][0] == 'I':
	140	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
	141
	142	def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
[939]	143	'Needs a doc string'
[926]	144	atxIds = ['dAx:','dAy:','dAz:']
	145	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
	146	TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
	147	LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
	148	SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
	149	PIds = ['Px:','Py:','Pz:']
	150	OIds = ['Oa:','Oi:','Oj:','Ok:']
	151	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
	152	if not RBIds['Vector'] and not RBIds['Residue']:
	153	return
	154	VRBIds = RBIds['Vector']
	155	RRBIds = RBIds['Residue']
	156	RBData = rigidbodyDict
	157	for item in parmDict:
	158	if 'RB' in item:
	159	dFdvDict[item] = 0. #NB: this is a vector which is no. refl. long & must be filled!
	160	General = Phase['General']
	161	cell = General['Cell'][1:7]
	162	Amat,Bmat = G2lat.cell2AB(cell)
	163	rpd = np.pi/180.
	164	rpd2 = rpd**2
	165	g = nl.inv(np.inner(Bmat,Bmat))
	166	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
	167	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
	168	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
	169	pfx = str(Phase['pId'])+'::'
	170	RBModels = Phase['RBModels']
	171	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
	172	VModel = RBData['Vector'][RBObj['RBId']]
	173	Q = RBObj['Orient'][0]
	174	Pos = RBObj['Orig'][0]
	175	jrb = VRBIds.index(RBObj['RBId'])
	176	rbsx = str(irb)+':'+str(jrb)
	177	dXdv = []
	178	for iv in range(len(VModel['VectMag'])):
	179	dCdv = []
	180	for vec in VModel['rbVect'][iv]:
	181	dCdv.append(G2mth.prodQVQ(Q,vec))
	182	dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
	183	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
	184	for ia,atId in enumerate(RBObj['Ids']):
	185	atNum = AtLookup[atId]
[953]	186	dx = 0.00001
[926]	187	for iv in range(len(VModel['VectMag'])):
	188	for ix in [0,1,2]:
	189	dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	190	for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
	191	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
[954]	192	for iv in range(4):
[926]	193	Q[iv] -= dx
[954]	194	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
[926]	195	Q[iv] += 2.*dx
[954]	196	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
[926]	197	Q[iv] -= dx
	198	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
	199	for ix in [0,1,2]:
	200	dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	201	X = G2mth.prodQVQ(Q,Cart[ia])
	202	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
	203	dFdu = G2lat.U6toUij(dFdu.T)
	204	dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
	205	dFdu = G2lat.UijtoU6(dFdu)
	206	atNum = AtLookup[atId]
	207	if 'T' in RBObj['ThermalMotion'][0]:
	208	for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
	209	dFdvDict[pfx+name+rbsx] += dFdu[i]
	210	if 'L' in RBObj['ThermalMotion'][0]:
	211	dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
	212	dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
	213	dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
	214	dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
	215	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
	216	dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2(-dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
	217	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
	218	dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2(-dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
	219	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
	220	if 'S' in RBObj['ThermalMotion'][0]:
	221	dFdvDict[pfx+'RBVS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
	222	dFdvDict[pfx+'RBVS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
	223	dFdvDict[pfx+'RBVS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
	224	dFdvDict[pfx+'RBVS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
	225	dFdvDict[pfx+'RBVS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
	226	dFdvDict[pfx+'RBVS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
	227	dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
	228	dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
	229	if 'U' in RBObj['ThermalMotion'][0]:
	230	dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
	231
	232
	233	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
	234	Q = RBObj['Orient'][0]
	235	Pos = RBObj['Orig'][0]
	236	jrb = RRBIds.index(RBObj['RBId'])
	237	torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
	238	rbsx = str(irb)+':'+str(jrb)
	239	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
	240	for itors,tors in enumerate(RBObj['Torsions']): #derivative error?
	241	tname = pfx+'RBRTr;'+str(itors)+':'+rbsx
	242	orId,pvId = torData[itors][:2]
	243	pivotVec = Cart[orId]-Cart[pvId]
	244	QA = G2mth.AVdeg2Q(-0.001,pivotVec)
	245	QB = G2mth.AVdeg2Q(0.001,pivotVec)
	246	for ir in torData[itors][3]:
	247	atNum = AtLookup[RBObj['Ids'][ir]]
	248	rVec = Cart[ir]-Cart[pvId]
	249	dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
	250	dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
	251	for ix in [0,1,2]:
	252	dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	253	for ia,atId in enumerate(RBObj['Ids']):
	254	atNum = AtLookup[atId]
[954]	255	dx = 0.00001
[926]	256	for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
	257	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
	258	for iv in range(4):
	259	Q[iv] -= dx
[954]	260	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
[926]	261	Q[iv] += 2.*dx
[954]	262	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
[926]	263	Q[iv] -= dx
	264	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
	265	for ix in [0,1,2]:
	266	dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
	267	X = G2mth.prodQVQ(Q,Cart[ia])
	268	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
	269	dFdu = G2lat.U6toUij(dFdu.T)
	270	dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
	271	dFdu = G2lat.UijtoU6(dFdu)
	272	atNum = AtLookup[atId]
	273	if 'T' in RBObj['ThermalMotion'][0]:
	274	for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
	275	dFdvDict[pfx+name+rbsx] += dFdu[i]
	276	if 'L' in RBObj['ThermalMotion'][0]:
	277	dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
	278	dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
	279	dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
	280	dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
	281	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
	282	dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2(dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
	283	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
	284	dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2(dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
	285	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
	286	if 'S' in RBObj['ThermalMotion'][0]:
	287	dFdvDict[pfx+'RBRS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
	288	dFdvDict[pfx+'RBRS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
	289	dFdvDict[pfx+'RBRS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
	290	dFdvDict[pfx+'RBRS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
	291	dFdvDict[pfx+'RBRS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
	292	dFdvDict[pfx+'RBRS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
	293	dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
	294	dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
	295	if 'U' in RBObj['ThermalMotion'][0]:
	296	dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
	297
	298	################################################################################
	299	##### Penalty & restraint functions
	300	################################################################################
	301
	302	def penaltyFxn(HistoPhases,parmDict,varyList):
[939]	303	'Needs a doc string'
[926]	304	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	305	pNames = []
	306	pVals = []
	307	pWt = []
	308	negWt = {}
	309	for phase in Phases:
	310	pId = Phases[phase]['pId']
	311	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
	312	General = Phases[phase]['General']
	313	textureData = General['SH Texture']
	314	SGData = General['SGData']
	315	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
	316	cell = General['Cell'][1:7]
	317	Amat,Bmat = G2lat.cell2AB(cell)
	318	if phase not in restraintDict:
	319	continue
	320	phaseRest = restraintDict[phase]
	321	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
	322	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
	323	['ChemComp','Sites'],['Texture','HKLs']]
	324	for name,rest in names:
	325	itemRest = phaseRest[name]
	326	if itemRest[rest] and itemRest['Use']:
	327	wt = itemRest['wtFactor']
	328	if name in ['Bond','Angle','Plane','Chiral']:
	329	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
	330	pNames.append(str(pId)+':'+name+':'+str(i))
	331	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	332	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
	333	if name == 'Bond':
	334	calc = G2mth.getRestDist(XYZ,Amat)
	335	elif name == 'Angle':
	336	calc = G2mth.getRestAngle(XYZ,Amat)
	337	elif name == 'Plane':
	338	calc = G2mth.getRestPlane(XYZ,Amat)
	339	elif name == 'Chiral':
	340	calc = G2mth.getRestChiral(XYZ,Amat)
	341	pVals.append(obs-calc)
	342	pWt.append(wt/esd**2)
	343	elif name in ['Torsion','Rama']:
	344	coeffDict = itemRest['Coeff']
	345	for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
	346	pNames.append(str(pId)+':'+name+':'+str(i))
	347	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	348	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
	349	if name == 'Torsion':
	350	tor = G2mth.getRestTorsion(XYZ,Amat)
	351	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
	352	else:
	353	phi,psi = G2mth.getRestRama(XYZ,Amat)
	354	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
	355	pVals.append(obs-calc)
	356	pWt.append(wt/esd**2)
	357	elif name == 'ChemComp':
	358	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
	359	pNames.append(str(pId)+':'+name+':'+str(i))
	360	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
	361	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
	362	calc = np.sum(mulfracfactors)
	363	pVals.append(obs-calc)
	364	pWt.append(wt/esd**2)
	365	elif name == 'Texture':
	366	SHkeys = textureData['SH Coeff'][1].keys()
	367	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
	368	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
	369	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
	370	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
	371	PH = np.array(hkl)
	372	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
	373	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
	374	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
	375	Z1 = -ma.masked_greater(Z,0.0)
	376	IndZ1 = np.array(ma.nonzero(Z1))
	377	for ind in IndZ1.T:
	378	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
	379	pVals.append(Z1[ind[0]][ind[1]])
	380	pWt.append(wt/esd1**2)
	381	if ifesd2:
	382	Z2 = 1.-Z
	383	for ind in np.ndindex(grid,grid):
	384	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
	385	pVals.append(Z1[ind[0]][ind[1]])
	386	pWt.append(wt/esd2**2)
	387
	388	for item in varyList:
	389	if 'PWLref' in item and parmDict[item] < 0.:
	390	pId = int(item.split(':')[0])
	391	if negWt[pId]:
	392	pNames.append(item)
	393	pVals.append(-parmDict[item])
	394	pWt.append(negWt[pId])
	395	pVals = np.array(pVals)
	396	pWt = np.array(pWt) #should this be np.sqrt?
	397	return pNames,pVals,pWt
	398
	399	def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
[939]	400	'Needs a doc string'
[926]	401	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	402	pDerv = np.zeros((len(varyList),len(pVal)))
	403	for phase in Phases:
	404	# if phase not in restraintDict:
	405	# continue
	406	pId = Phases[phase]['pId']
	407	General = Phases[phase]['General']
	408	SGData = General['SGData']
	409	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
	410	cell = General['Cell'][1:7]
	411	Amat,Bmat = G2lat.cell2AB(cell)
	412	textureData = General['SH Texture']
	413
	414	SHkeys = textureData['SH Coeff'][1].keys()
	415	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
	416	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
	417	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
	418	sam = SamSym[textureData['Model']]
	419	phaseRest = restraintDict.get(phase,{})
	420	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
	421	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
	422	'ChemComp':'Sites','Texture':'HKLs'}
	423	lasthkl = np.array([0,0,0])
	424	for ip,pName in enumerate(pNames):
	425	pnames = pName.split(':')
	426	if pId == int(pnames[0]):
	427	name = pnames[1]
	428	if 'PWL' in pName:
	429	pDerv[varyList.index(pName)][ip] += 1.
	430	continue
	431	id = int(pnames[2])
	432	itemRest = phaseRest[name]
	433	if name in ['Bond','Angle','Plane','Chiral']:
	434	indx,ops,obs,esd = itemRest[names[name]][id]
	435	dNames = []
	436	for ind in indx:
	437	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
	438	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	439	if name == 'Bond':
	440	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
	441	elif name == 'Angle':
	442	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
	443	elif name == 'Plane':
	444	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
	445	elif name == 'Chiral':
	446	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
	447	elif name in ['Torsion','Rama']:
	448	coffDict = itemRest['Coeff']
	449	indx,ops,cofName,esd = itemRest[names[name]][id]
	450	dNames = []
	451	for ind in indx:
	452	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
	453	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
	454	if name == 'Torsion':
	455	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
	456	else:
	457	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
	458	elif name == 'ChemComp':
	459	indx,factors,obs,esd = itemRest[names[name]][id]
	460	dNames = []
	461	for ind in indx:
	462	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
	463	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
	464	deriv = mul*factors
	465	elif 'Texture' in name:
	466	deriv = []
	467	dNames = []
	468	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
	469	hkl = np.array(hkl)
	470	if np.any(lasthkl-hkl):
	471	PH = np.array(hkl)
	472	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
	473	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
	474	lasthkl = copy.copy(hkl)
	475	if 'unit' in name:
	476	pass
	477	else:
	478	gam = float(pnames[3])
	479	psi = float(pnames[4])
	480	for SHname in ODFln:
	481	l,m,n = eval(SHname[1:])
	482	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
	483	dNames += [str(pId)+'::'+SHname]
	484	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
	485	for dName,drv in zip(dNames,deriv):
	486	try:
	487	ind = varyList.index(dName)
	488	pDerv[ind][ip] += drv
	489	except ValueError:
	490	pass
	491	return pDerv
	492
	493	################################################################################
	494	##### Function & derivative calculations
	495	################################################################################
	496
	497	def GetAtomFXU(pfx,calcControls,parmDict):
[939]	498	'Needs a doc string'
[926]	499	Natoms = calcControls['Natoms'][pfx]
	500	Tdata = Natoms*[' ',]
	501	Mdata = np.zeros(Natoms)
	502	IAdata = Natoms*[' ',]
	503	Fdata = np.zeros(Natoms)
	504	FFdata = []
	505	BLdata = []
	506	Xdata = np.zeros((3,Natoms))
	507	dXdata = np.zeros((3,Natoms))
	508	Uisodata = np.zeros(Natoms)
	509	Uijdata = np.zeros((6,Natoms))
	510	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
	511	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
	512	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
	513	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
	514	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
	515	for iatm in range(Natoms):
	516	for key in keys:
	517	parm = pfx+key+str(iatm)
	518	if parm in parmDict:
	519	keys[key][iatm] = parmDict[parm]
	520	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
	521
	522	def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
	523	''' Compute structure factors for all h,k,l for phase
[939]	524	puts the result, F^2, in each ref[8] in refList
[926]	525	input:
[939]	526
	527	:param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
	528	:param np.array G: reciprocal metric tensor
	529	:param str pfx: phase id string
	530	:param dict SGData: space group info. dictionary output from SpcGroup
	531	:param dict calcControls:
	532	:param dict ParmDict:
	533
[926]	534	'''
	535	twopi = 2.0*np.pi
	536	twopisq = 2.0np.pi*2
	537	phfx = pfx.split(':')[0]+hfx
	538	ast = np.sqrt(np.diag(G))
	539	Mast = twopisq*np.multiply.outer(ast,ast)
	540	FFtables = calcControls['FFtables']
	541	BLtables = calcControls['BLtables']
	542	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
	543	FF = np.zeros(len(Tdata))
[1078]	544	if 'N' in calcControls[hfx+'histType']:
[926]	545	FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
	546	else:
	547	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	548	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	549	Uij = np.array(G2lat.U6toUij(Uijdata))
	550	bij = Mast*Uij.T
	551	for refl in refList:
	552	fbs = np.array([0,0])
	553	H = refl[:3]
	554	SQ = 1./(2.refl[4])*2
	555	SQfactor = 4.0SQtwopisq
	556	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
	557	if not len(refl[-1]): #no form factors
[1078]	558	if 'N' in calcControls[hfx+'histType']:
[942]	559	refl[-1] = G2el.getBLvalues(BLtables)
[926]	560	else: #'X'
[942]	561	refl[-1] = G2el.getFFvalues(FFtables,SQ)
[926]	562	for i,El in enumerate(Tdata):
	563	FF[i] = refl[-1][El]
	564	Uniq = refl[11]
	565	phi = refl[12]
	566	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
	567	sinp = np.sin(phase)
	568	cosp = np.cos(phase)
	569	occ = Mdata*Fdata/len(Uniq)
	570	biso = -SQfactor*Uisodata
	571	Tiso = np.where(biso<1.,np.exp(biso),1.0)
	572	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
	573	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
	574	Tcorr = Tiso*Tuij
	575	fa = np.array([(FF+FP-Bab)occcospTcorr,-FPPoccsinpTcorr])
	576	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
	577	if not SGData['SGInv']:
	578	fb = np.array([(FF+FP-Bab)occsinpTcorr,FPPocccospTcorr])
	579	fbs = np.sum(np.sum(fb,axis=1),axis=1)
	580	fasq = fas**2
	581	fbsq = fbs**2 #imaginary
	582	refl[9] = np.sum(fasq)+np.sum(fbsq)
	583	refl[10] = atan2d(fbs[0],fas[0])
	584
	585	def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
[939]	586	'Needs a doc string'
[926]	587	twopi = 2.0*np.pi
	588	twopisq = 2.0np.pi*2
	589	phfx = pfx.split(':')[0]+hfx
	590	ast = np.sqrt(np.diag(G))
	591	Mast = twopisq*np.multiply.outer(ast,ast)
	592	FFtables = calcControls['FFtables']
	593	BLtables = calcControls['BLtables']
	594	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
	595	FF = np.zeros(len(Tdata))
[1078]	596	if 'N' in calcControls[hfx+'histType']:
[926]	597	FP = 0.
	598	FPP = 0.
	599	else:
	600	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
	601	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
	602	Uij = np.array(G2lat.U6toUij(Uijdata))
	603	bij = Mast*Uij.T
	604	dFdvDict = {}
	605	dFdfr = np.zeros((len(refList),len(Mdata)))
	606	dFdx = np.zeros((len(refList),len(Mdata),3))
	607	dFdui = np.zeros((len(refList),len(Mdata)))
	608	dFdua = np.zeros((len(refList),len(Mdata),6))
	609	dFdbab = np.zeros((len(refList),2))
	610	for iref,refl in enumerate(refList):
	611	H = np.array(refl[:3])
	612	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
	613	SQfactor = 8.0SQnp.pi**2
	614	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
	615	Bab = parmDict[phfx+'BabA']*dBabdA
	616	for i,El in enumerate(Tdata):
	617	FF[i] = refl[-1][El]
	618	Uniq = refl[11]
	619	phi = refl[12]
	620	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
	621	sinp = np.sin(phase)
	622	cosp = np.cos(phase)
	623	occ = Mdata*Fdata/len(Uniq)
	624	biso = -SQfactor*Uisodata
	625	Tiso = np.where(biso<1.,np.exp(biso),1.0)
	626	HbH = -np.inner(H,np.inner(bij,H))
	627	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
	628	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
	629	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
	630	Tcorr = Tiso*Tuij
	631	fot = (FF+FP-Bab)occTcorr
	632	fotp = FPPoccTcorr
	633	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
	634	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
	635
	636	fas = np.sum(np.sum(fa,axis=1),axis=1)
	637	fbs = np.sum(np.sum(fb,axis=1),axis=1)
	638	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
	639	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
	640	#sum below is over Uniq
	641	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
	642	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
	643	dfadui = np.sum(-SQfactor*fa,axis=2)
	644	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
	645	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
	646	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
	647	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
	648	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
	649	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
	650	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
	651	dFdbab[iref] = np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
	652	if not SGData['SGInv']:
	653	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom
	654	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
	655	dfbdui = np.sum(-SQfactor*fb,axis=2)
	656	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
	657	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
	658	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
	659	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
	660	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
	661	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
	662	dFdbab[iref] += np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
	663	#loop over atoms - each dict entry is list of derivatives for all the reflections
	664	for i in range(len(Mdata)):
	665	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
	666	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
	667	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
	668	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
	669	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
	670	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
	671	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
	672	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
	673	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
	674	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
	675	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
	676	dFdvDict[pfx+'BabA'] = dFdbab.T[0]
	677	dFdvDict[pfx+'BabU'] = dFdbab.T[1]
	678	return dFdvDict
	679
	680	def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
	681	''' Single crystal extinction function; puts correction in ref[13] and returns
	682	corrections needed for derivatives
	683	'''
	684	ref[13] = 1.0
	685	dervCor = 1.0
	686	dervDict = {}
	687	if calcControls[phfx+'EType'] != 'None':
	688	cos2T = 1.0-0.5(parmDict[hfx+'Lam']/ref[4])*2 #cos(2theta)
	689	if 'SXC' in parmDict[hfx+'Type']:
	690	AV = 7.9406e5/parmDict[pfx+'Vol']**2
	691	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
	692	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
	693	elif 'SNT' in parmDict[hfx+'Type']:
	694	AV = 1.e7/parmDict[pfx+'Vol']**2
	695	PL = 1./(4.refl[4]*2)
	696	P12 = 1.0
	697	elif 'SNC' in parmDict[hfx+'Type']:
	698	AV = 1.e7/parmDict[pfx+'Vol']**2
	699	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
	700	P12 = 1.0
	701
	702	PLZ = AVP12parmDict[hfx+'Lam']*2ref[7]
	703	if 'Primary' in calcControls[phfx+'EType']:
	704	PLZ *= 1.5
	705	else:
	706	PLZ *= calcControls[phfx+'Tbar']
	707
	708	if 'Primary' in calcControls[phfx+'EType']:
	709	PSIG = parmDict[phfx+'Ep']
	710	elif 'I & II' in calcControls[phfx+'EType']:
	711	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
	712	elif 'Type II' in calcControls[phfx+'EType']:
	713	PSIG = parmDict[phfx+'Es']
	714	else: # 'Secondary Type I'
	715	PSIG = parmDict[phfx+'Eg']/PL
	716
	717	AG = 0.58+0.48cos2T+0.24cos2T**2
	718	AL = 0.025+0.285*cos2T
	719	BG = 0.02-0.025*cos2T
	720	BL = 0.15-0.2(0.75-cos2T)*2
	721	if cos2T < 0.:
	722	BL = -0.45*cos2T
	723	CG = 2.
	724	CL = 2.
	725	PF = PLZ*PSIG
	726
	727	if 'Gaussian' in calcControls[phfx+'EApprox']:
	728	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
	729	extCor = np.sqrt(PF4)
	730	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
	731	else:
	732	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
	733	extCor = np.sqrt(PF4)
	734	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
	735
	736	dervCor = (1.+PF)*PF3
	737	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
	738	dervDict[phfx+'Ep'] = -ref[7]PLZPF3
	739	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
	740	dervDict[phfx+'Es'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
	741	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
	742	dervDict[phfx+'Eg'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
	743
	744	ref[13] = 1./extCor
	745	return dervCor,dervDict
	746
	747
	748	def Dict2Values(parmdict, varylist):
	749	'''Use before call to leastsq to setup list of values for the parameters
	750	in parmdict, as selected by key in varylist'''
	751	return [parmdict[key] for key in varylist]
	752
	753	def Values2Dict(parmdict, varylist, values):
	754	''' Use after call to leastsq to update the parameter dictionary with
	755	values corresponding to keys in varylist'''
	756	parmdict.update(zip(varylist,values))
	757
	758	def GetNewCellParms(parmDict,varyList):
[939]	759	'Needs a doc string'
[926]	760	newCellDict = {}
	761	Anames = ['A'+str(i) for i in range(6)]
	762	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
	763	for item in varyList:
	764	keys = item.split(':')
	765	if keys[2] in Ddict:
	766	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
	767	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
	768	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
	769	return newCellDict # is e.g. {'0::D11':A0+D11}
	770
	771	def ApplyXYZshifts(parmDict,varyList):
	772	'''
[939]	773	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
	774
	775	:param dict parmDict: parameter dictionary
	776	:param list varyList: list of variables
	777	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
	778
	779	'''
[926]	780	newAtomDict = {}
	781	for item in parmDict:
	782	if 'dA' in item:
	783	parm = ''.join(item.split('d'))
	784	parmDict[parm] += parmDict[item]
	785	newAtomDict[item] = [parm,parmDict[parm]]
	786	return newAtomDict
	787
	788	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
[939]	789	'Spherical harmonics texture'
[926]	790	IFCoup = 'Bragg' in calcControls[hfx+'instType']
	791	odfCor = 1.0
	792	H = refl[:3]
	793	cell = G2lat.Gmat2cell(g)
	794	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
	795	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
	796	phi,beta = G2lat.CrsAng(H,cell,SGData)
	797	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
	798	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
	799	for item in SHnames:
	800	L,M,N = eval(item.strip('C'))
	801	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
	802	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
	803	Lnorm = G2lat.Lnorm(L)
	804	odfCor += parmDict[pfx+item]LnormKcl*Ksl
	805	return odfCor
	806
	807	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
[939]	808	'Spherical harmonics texture derivatives'
[926]	809	FORPI = 4.0*np.pi
	810	IFCoup = 'Bragg' in calcControls[hfx+'instType']
	811	odfCor = 1.0
	812	dFdODF = {}
	813	dFdSA = [0,0,0]
	814	H = refl[:3]
	815	cell = G2lat.Gmat2cell(g)
	816	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
	817	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
	818	phi,beta = G2lat.CrsAng(H,cell,SGData)
	819	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
	820	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
	821	for item in SHnames:
	822	L,M,N = eval(item.strip('C'))
	823	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
	824	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
	825	Lnorm = G2lat.Lnorm(L)
	826	odfCor += parmDict[pfx+item]LnormKcl*Ksl
	827	dFdODF[pfx+item] = LnormKclKsl
	828	for i in range(3):
	829	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
	830	return odfCor,dFdODF,dFdSA
	831
	832	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
[939]	833	'spherical harmonics preferred orientation (cylindrical symmetry only)'
[926]	834	odfCor = 1.0
	835	H = refl[:3]
	836	cell = G2lat.Gmat2cell(g)
	837	Sangl = [0.,0.,0.]
	838	if 'Bragg' in calcControls[hfx+'instType']:
	839	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
	840	IFCoup = True
	841	else:
	842	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
	843	IFCoup = False
	844	phi,beta = G2lat.CrsAng(H,cell,SGData)
	845	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
	846	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
	847	for item in SHnames:
	848	L,N = eval(item.strip('C'))
[1039]	849	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
[926]	850	odfCor += parmDict[phfx+item]LnormKcsl
[1039]	851	return np.squeeze(odfCor)
[926]	852
	853	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
[939]	854	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
[926]	855	FORPI = 12.5663706143592
	856	odfCor = 1.0
	857	dFdODF = {}
	858	H = refl[:3]
	859	cell = G2lat.Gmat2cell(g)
	860	Sangl = [0.,0.,0.]
	861	if 'Bragg' in calcControls[hfx+'instType']:
	862	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
	863	IFCoup = True
	864	else:
	865	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
	866	IFCoup = False
	867	phi,beta = G2lat.CrsAng(H,cell,SGData)
	868	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
	869	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
	870	for item in SHnames:
	871	L,N = eval(item.strip('C'))
	872	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
	873	odfCor += parmDict[phfx+item]LnormKcsl
	874	dFdODF[phfx+item] = Kcsl*Lnorm
	875	return odfCor,dFdODF
	876
	877	def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
[939]	878	'Needs a doc string'
[926]	879	POcorr = 1.0
	880	if calcControls[phfx+'poType'] == 'MD':
	881	MD = parmDict[phfx+'MD']
	882	if MD != 1.0:
	883	MDAxis = calcControls[phfx+'MDAxis']
	884	sumMD = 0
	885	for H in refl[11]:
	886	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
	887	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
	888	sumMD += A**3
	889	POcorr = sumMD/len(refl[11])
	890	else: #spherical harmonics
	891	if calcControls[phfx+'SHord']:
	892	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
	893	return POcorr
	894
	895	def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
[939]	896	'Needs a doc string'
[926]	897	POcorr = 1.0
	898	POderv = {}
	899	if calcControls[phfx+'poType'] == 'MD':
	900	MD = parmDict[phfx+'MD']
	901	MDAxis = calcControls[phfx+'MDAxis']
	902	sumMD = 0
	903	sumdMD = 0
	904	for H in refl[11]:
	905	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
	906	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
	907	sumMD += A**3
	908	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
	909	POcorr = sumMD/len(refl[11])
	910	POderv[phfx+'MD'] = sumdMD/len(refl[11])
	911	else: #spherical harmonics
	912	if calcControls[phfx+'SHord']:
	913	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
	914	return POcorr,POderv
	915
	916	def GetAbsorb(refl,hfx,calcControls,parmDict):
[939]	917	'Needs a doc string'
[926]	918	if 'Debye' in calcControls[hfx+'instType']:
	919	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
	920	else:
	921	return 1.0
	922
	923	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
[939]	924	'Needs a doc string'
[926]	925	if 'Debye' in calcControls[hfx+'instType']:
	926	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
	927	else:
	928	return 0.0
	929
	930	def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
[939]	931	'Needs a doc string'
[926]	932	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
	933	if 'X' in parmDict[hfx+'Type']:
	934	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
	935	Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
	936	if pfx+'SHorder' in parmDict:
	937	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
	938	Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
	939	refl[13] = Icorr
	940
	941	def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
[939]	942	'Needs a doc string'
[926]	943	dIdsh = 1./parmDict[hfx+'Scale']
	944	dIdsp = 1./parmDict[phfx+'Scale']
	945	if 'X' in parmDict[hfx+'Type']:
	946	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
	947	dIdPola /= pola
	948	else: #'N'
	949	dIdPola = 0.0
	950	POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
	951	for iPO in dIdPO:
	952	dIdPO[iPO] /= POcorr
	953	dFdODF = {}
	954	dFdSA = [0,0,0]
	955	if pfx+'SHorder' in parmDict:
	956	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
	957	for iSH in dFdODF:
	958	dFdODF[iSH] /= odfCor
	959	for i in range(3):
	960	dFdSA[i] /= odfCor
	961	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
	962	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
	963
	964	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
[939]	965	'Needs a doc string'
[926]	966	costh = cosd(refl[5]/2.)
	967	#crystallite size
	968	if calcControls[phfx+'SizeType'] == 'isotropic':
	969	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
	970	elif calcControls[phfx+'SizeType'] == 'uniaxial':
	971	H = np.array(refl[:3])
	972	P = np.array(calcControls[phfx+'SizeAxis'])
	973	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	974	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
	975	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
	976	else: #ellipsoidal crystallites
	977	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
	978	H = np.array(refl[:3])
	979	lenR = G2pwd.ellipseSize(H,Sij,GB)
	980	Sgam = 1.8wave/(np.picosth*lenR)
	981	#microstrain
	982	if calcControls[phfx+'MustrainType'] == 'isotropic':
	983	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
	984	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
	985	H = np.array(refl[:3])
	986	P = np.array(calcControls[phfx+'MustrainAxis'])
	987	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	988	Si = parmDict[phfx+'Mustrain;i']
	989	Sa = parmDict[phfx+'Mustrain;a']
	990	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
	991	else: #generalized - P.W. Stephens model
	992	pwrs = calcControls[phfx+'MuPwrs']
	993	sum = 0
	994	for i,pwr in enumerate(pwrs):
	995	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
	996	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
	997	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
	998	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
	999	sig /= ateln2
	1000	return sig,gam
	1001
	1002	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
[939]	1003	'Needs a doc string'
[926]	1004	gamDict = {}
	1005	sigDict = {}
	1006	costh = cosd(refl[5]/2.)
	1007	tanth = tand(refl[5]/2.)
	1008	#crystallite size derivatives
	1009	if calcControls[phfx+'SizeType'] == 'isotropic':
	1010	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
	1011	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
	1012	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
	1013	elif calcControls[phfx+'SizeType'] == 'uniaxial':
	1014	H = np.array(refl[:3])
	1015	P = np.array(calcControls[phfx+'SizeAxis'])
	1016	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1017	Si = parmDict[phfx+'Size;i']
	1018	Sa = parmDict[phfx+'Size;a']
	1019	gami = (1.8wave/np.pi)/(SiSa)
	1020	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
	1021	Sgam = gami*sqtrm
	1022	gam = Sgam/costh
	1023	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
	1024	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
	1025	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
	1026	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
	1027	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1028	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1029	else: #ellipsoidal crystallites
	1030	const = 1.8wave/(np.picosth)
	1031	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
	1032	H = np.array(refl[:3])
	1033	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
	1034	Sgam = 1.8wave/(np.picosth*lenR)
	1035	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
	1036	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
	1037	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
	1038	gamDict[phfx+'Size;mx'] = Sgam
	1039	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
	1040
	1041	#microstrain derivatives
	1042	if calcControls[phfx+'MustrainType'] == 'isotropic':
	1043	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
	1044	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
	1045	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
	1046	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
	1047	H = np.array(refl[:3])
	1048	P = np.array(calcControls[phfx+'MustrainAxis'])
	1049	cosP,sinP = G2lat.CosSinAngle(H,P,G)
	1050	Si = parmDict[phfx+'Mustrain;i']
	1051	Sa = parmDict[phfx+'Mustrain;a']
	1052	gami = 0.018SiSa*tanth/np.pi
	1053	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
	1054	Mgam = gami/sqtrm
	1055	dsi = -gamiSicosP2/sqtrm3
	1056	dsa = -gamiSasinP2/sqtrm3
	1057	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
	1058	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
	1059	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1060	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1061	else: #generalized - P.W. Stephens model
	1062	pwrs = calcControls[phfx+'MuPwrs']
	1063	const = 0.018refl[4]2tanth
	1064	sum = 0
	1065	for i,pwr in enumerate(pwrs):
	1066	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
	1067	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
	1068	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
	1069	sigDict[phfx+'Mustrain:'+str(i)] = \
	1070	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
	1071	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
	1072	for i in range(len(pwrs)):
	1073	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
	1074	gamDict[phfx+'Mustrain;mx'] = Mgam
	1075	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
	1076	return sigDict,gamDict
	1077
	1078	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
[939]	1079	'Needs a doc string'
[926]	1080	h,k,l = refl[:3]
	1081	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
	1082	d = np.sqrt(dsq)
	1083
	1084	refl[4] = d
	1085	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
	1086	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
	1087	if 'Bragg' in calcControls[hfx+'instType']:
	1088	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
	1089	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
	1090	else: #Debye-Scherrer - simple but maybe not right
	1091	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
	1092	return pos
	1093
	1094	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
[939]	1095	'Needs a doc string'
[926]	1096	dpr = 180./np.pi
	1097	h,k,l = refl[:3]
	1098	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
	1099	dst = np.sqrt(dstsq)
	1100	pos = refl[5]-parmDict[hfx+'Zero']
	1101	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
	1102	dpdw = const*dst
	1103	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
	1104	dpdA = constwave/(2.0*dst)
	1105	dpdZ = 1.0
	1106	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
	1107	if 'Bragg' in calcControls[hfx+'instType']:
	1108	dpdSh = -4.constcosd(pos/2.0)
	1109	dpdTr = -constsind(pos)100.0
	1110	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
	1111	else: #Debye-Scherrer - simple but maybe not right
	1112	dpdXd = -const*cosd(pos)
	1113	dpdYd = -const*sind(pos)
	1114	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
	1115
	1116	def GetHStrainShift(refl,SGData,phfx,parmDict):
[939]	1117	'Needs a doc string'
[926]	1118	laue = SGData['SGLaue']
	1119	uniq = SGData['SGUniq']
	1120	h,k,l = refl[:3]
	1121	if laue in ['m3','m3m']:
	1122	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
	1123	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
	1124	elif laue in ['6/m','6/mmm','3m1','31m','3']:
	1125	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
	1126	elif laue in ['3R','3mR']:
	1127	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
	1128	elif laue in ['4/m','4/mmm']:
	1129	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
	1130	elif laue in ['mmm']:
	1131	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
	1132	elif laue in ['2/m']:
	1133	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
	1134	if uniq == 'a':
	1135	Dij += parmDict[phfx+'D23']kl
	1136	elif uniq == 'b':
	1137	Dij += parmDict[phfx+'D13']hl
	1138	elif uniq == 'c':
	1139	Dij += parmDict[phfx+'D12']hk
	1140	else:
	1141	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
	1142	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
	1143	return -Dijrefl[4]2tand(refl[5]/2.0)
	1144
	1145	def GetHStrainShiftDerv(refl,SGData,phfx):
[939]	1146	'Needs a doc string'
[926]	1147	laue = SGData['SGLaue']
	1148	uniq = SGData['SGUniq']
	1149	h,k,l = refl[:3]
	1150	if laue in ['m3','m3m']:
	1151	dDijDict = {phfx+'D11':h2+k2+l**2,
	1152	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
	1153	elif laue in ['6/m','6/mmm','3m1','31m','3']:
	1154	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
	1155	elif laue in ['3R','3mR']:
	1156	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
	1157	elif laue in ['4/m','4/mmm']:
	1158	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
	1159	elif laue in ['mmm']:
	1160	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
	1161	elif laue in ['2/m']:
	1162	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
	1163	if uniq == 'a':
	1164	dDijDict[phfx+'D23'] = k*l
	1165	elif uniq == 'b':
	1166	dDijDict[phfx+'D13'] = h*l
	1167	elif uniq == 'c':
	1168	dDijDict[phfx+'D12'] = h*k
	1169	names.append()
	1170	else:
	1171	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
	1172	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
	1173	for item in dDijDict:
	1174	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
	1175	return dDijDict
	1176
	1177	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
[939]	1178	'Needs a doc string'
[926]	1179	histoList = Histograms.keys()
	1180	histoList.sort()
	1181	for histogram in histoList:
	1182	if 'PWDR' in histogram[:4]:
	1183	Histogram = Histograms[histogram]
	1184	hId = Histogram['hId']
	1185	hfx = ':%d:'%(hId)
	1186	Limits = calcControls[hfx+'Limits']
	1187	shl = max(parmDict[hfx+'SH/L'],0.0005)
	1188	Ka2 = False
	1189	kRatio = 0.0
	1190	if hfx+'Lam1' in parmDict.keys():
	1191	Ka2 = True
	1192	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
	1193	kRatio = parmDict[hfx+'I(L2)/I(L1)']
	1194	x,y,w,yc,yb,yd = Histogram['Data']
	1195	xB = np.searchsorted(x,Limits[0])
	1196	xF = np.searchsorted(x,Limits[1])
	1197	ymb = np.array(y-yb)
	1198	ymb = np.where(ymb,ymb,1.0)
	1199	ycmb = np.array(yc-yb)
	1200	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
	1201	refLists = Histogram['Reflection Lists']
	1202	for phase in refLists:
	1203	Phase = Phases[phase]
	1204	pId = Phase['pId']
	1205	phfx = '%d:%d:'%(pId,hId)
	1206	refList = refLists[phase]
	1207	sumFo = 0.0
	1208	sumdF = 0.0
	1209	sumFosq = 0.0
	1210	sumdFsq = 0.0
	1211	for refl in refList:
	1212	if 'C' in calcControls[hfx+'histType']:
	1213	yp = np.zeros_like(yb)
	1214	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
	1215	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
	1216	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
	1217	iFin2 = iFin
	1218	yp[iBeg:iFin] = refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
	1219	if Ka2:
	1220	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
	1221	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
	1222	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
	1223	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
	1224	if not iBeg2+iFin2: #peak below low limit - skip peak
	1225	continue
	1226	elif not iBeg2-iFin2: #peak above high limit - done
	1227	break
	1228	yp[iBeg2:iFin2] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
	1229	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[13](1.+kRatio)),0.0))
	1230	elif 'T' in calcControls[hfx+'histType']:
	1231	print 'TOF Undefined at present'
	1232	raise Exception #no TOF yet
	1233	Fo = np.sqrt(np.abs(refl[8]))
	1234	Fc = np.sqrt(np.abs(refl[9]))
	1235	sumFo += Fo
	1236	sumFosq += refl[8]**2
	1237	sumdF += np.abs(Fo-Fc)
	1238	sumdFsq += (refl[8]-refl[9])**2
[1030]	1239	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
	1240	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
	1241	Histogram['Residuals'][phfx+'Nref'] = len(refList)
[1035]	1242	Histogram['Residuals']['hId'] = hId
[1036]	1243	elif 'HKLF' in histogram[:4]:
[1038]	1244	Histogram = Histograms[histogram]
[1036]	1245	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
[926]	1246
	1247	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
[939]	1248	'Needs a doc string'
[926]	1249
	1250	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
	1251	U = parmDict[hfx+'U']
	1252	V = parmDict[hfx+'V']
	1253	W = parmDict[hfx+'W']
	1254	X = parmDict[hfx+'X']
	1255	Y = parmDict[hfx+'Y']
	1256	tanPos = tand(refl[5]/2.0)
	1257	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
	1258	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
	1259	sig = max(0.001,sig)
	1260	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
	1261	gam = max(0.001,gam)
	1262	return sig,gam
	1263
	1264	hId = Histogram['hId']
	1265	hfx = ':%d:'%(hId)
	1266	bakType = calcControls[hfx+'bakType']
	1267	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
	1268	yc = np.zeros_like(yb)
	1269
	1270	if 'C' in calcControls[hfx+'histType']:
	1271	shl = max(parmDict[hfx+'SH/L'],0.002)
	1272	Ka2 = False
	1273	if hfx+'Lam1' in parmDict.keys():
	1274	wave = parmDict[hfx+'Lam1']
	1275	Ka2 = True
	1276	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
	1277	kRatio = parmDict[hfx+'I(L2)/I(L1)']
	1278	else:
	1279	wave = parmDict[hfx+'Lam']
	1280	else:
	1281	print 'TOF Undefined at present'
	1282	raise ValueError
	1283	for phase in Histogram['Reflection Lists']:
	1284	refList = Histogram['Reflection Lists'][phase]
	1285	Phase = Phases[phase]
	1286	pId = Phase['pId']
	1287	pfx = '%d::'%(pId)
	1288	phfx = '%d:%d:'%(pId,hId)
	1289	hfx = ':%d:'%(hId)
	1290	SGData = Phase['General']['SGData']
	1291	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	1292	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	1293	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
	1294	Vst = np.sqrt(nl.det(G)) #V*
	1295	if not Phase['General'].get('doPawley'):
	1296	time0 = time.time()
	1297	StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
	1298	# print 'sf calc time: %.3fs'%(time.time()-time0)
	1299	time0 = time.time()
	1300	for refl in refList:
	1301	if 'C' in calcControls[hfx+'histType']:
	1302	h,k,l = refl[:3]
	1303	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
	1304	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
	1305	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
	1306	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
	1307	GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]
	1308	refl[13] = VstLorenz
	1309	if Phase['General'].get('doPawley'):
	1310	try:
	1311	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
	1312	refl[9] = parmDict[pInd]
	1313	except KeyError:
	1314	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
	1315	continue
	1316	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
	1317	iBeg = np.searchsorted(x,refl[5]-fmin)
	1318	iFin = np.searchsorted(x,refl[5]+fmax)
	1319	if not iBeg+iFin: #peak below low limit - skip peak
	1320	continue
	1321	elif not iBeg-iFin: #peak above high limit - done
	1322	break
[1025]	1323	yc[iBeg:iFin] += refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
[926]	1324	if Ka2:
	1325	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
	1326	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
	1327	iBeg = np.searchsorted(x,pos2-fmin)
	1328	iFin = np.searchsorted(x,pos2+fmax)
	1329	if not iBeg+iFin: #peak below low limit - skip peak
	1330	continue
	1331	elif not iBeg-iFin: #peak above high limit - done
	1332	return yc,yb
[1025]	1333	yc[iBeg:iFin] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here
[926]	1334	elif 'T' in calcControls[hfx+'histType']:
	1335	print 'TOF Undefined at present'
	1336	raise Exception #no TOF yet
	1337	# print 'profile calc time: %.3fs'%(time.time()-time0)
	1338	return yc,yb
	1339
	1340	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
[939]	1341	'Needs a doc string'
[926]	1342
	1343	def cellVaryDerv(pfx,SGData,dpdA):
	1344	if SGData['SGLaue'] in ['-1',]:
	1345	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
	1346	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
	1347	elif SGData['SGLaue'] in ['2/m',]:
	1348	if SGData['SGUniq'] == 'a':
	1349	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
	1350	elif SGData['SGUniq'] == 'b':
	1351	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
	1352	else:
	1353	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
	1354	elif SGData['SGLaue'] in ['mmm',]:
	1355	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
	1356	elif SGData['SGLaue'] in ['4/m','4/mmm']:
	1357	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
	1358	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
	1359	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
	1360	elif SGData['SGLaue'] in ['3R', '3mR']:
	1361	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
	1362	elif SGData['SGLaue'] in ['m3m','m3']:
	1363	return [[pfx+'A0',dpdA[0]]]
	1364
	1365	# create a list of dependent variables and set up a dictionary to hold their derivatives
	1366	dependentVars = G2mv.GetDependentVars()
	1367	depDerivDict = {}
	1368	for j in dependentVars:
	1369	depDerivDict[j] = np.zeros(shape=(len(x)))
	1370	#print 'dependent vars',dependentVars
	1371	lenX = len(x)
	1372	hId = Histogram['hId']
	1373	hfx = ':%d:'%(hId)
	1374	bakType = calcControls[hfx+'bakType']
	1375	dMdv = np.zeros(shape=(len(varylist),len(x)))
	1376	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
	1377	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
	1378	bBpos =varylist.index(hfx+'Back:0')
	1379	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
	1380	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
	1381	for name in varylist:
	1382	if 'Debye' in name:
	1383	id = int(name.split(':')[-1])
	1384	parm = name[:int(name.rindex(':'))]
	1385	ip = names.index(parm)
	1386	dMdv[varylist.index(name)] = dMddb[3*id+ip]
	1387	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
	1388	for name in varylist:
	1389	if 'BkPk' in name:
	1390	id = int(name.split(':')[-1])
	1391	parm = name[:int(name.rindex(':'))]
	1392	ip = names.index(parm)
	1393	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
	1394	cw = np.diff(x)
	1395	cw = np.append(cw,cw[-1])
	1396	if 'C' in calcControls[hfx+'histType']:
	1397	shl = max(parmDict[hfx+'SH/L'],0.002)
	1398	Ka2 = False
	1399	if hfx+'Lam1' in parmDict.keys():
	1400	wave = parmDict[hfx+'Lam1']
	1401	Ka2 = True
	1402	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
	1403	kRatio = parmDict[hfx+'I(L2)/I(L1)']
	1404	else:
	1405	wave = parmDict[hfx+'Lam']
	1406	else:
	1407	print 'TOF Undefined at present'
	1408	raise ValueError
	1409	for phase in Histogram['Reflection Lists']:
	1410	refList = Histogram['Reflection Lists'][phase]
	1411	Phase = Phases[phase]
	1412	SGData = Phase['General']['SGData']
	1413	pId = Phase['pId']
	1414	pfx = '%d::'%(pId)
	1415	phfx = '%d:%d:'%(pId,hId)
	1416	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	1417	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	1418	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
	1419	if not Phase['General'].get('doPawley'):
	1420	time0 = time.time()
	1421	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
	1422	# print 'sf-derv time %.3fs'%(time.time()-time0)
	1423	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
	1424	time0 = time.time()
	1425	for iref,refl in enumerate(refList):
	1426	if 'C' in calcControls[hfx+'histType']: #CW powder
	1427	h,k,l = refl[:3]
	1428	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
	1429	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
	1430	iBeg = np.searchsorted(x,refl[5]-fmin)
	1431	iFin = np.searchsorted(x,refl[5]+fmax)
	1432	if not iBeg+iFin: #peak below low limit - skip peak
	1433	continue
	1434	elif not iBeg-iFin: #peak above high limit - done
	1435	break
	1436	pos = refl[5]
	1437	tanth = tand(pos/2.0)
	1438	costh = cosd(pos/2.0)
	1439	lenBF = iFin-iBeg
	1440	dMdpk = np.zeros(shape=(6,lenBF))
[1025]	1441	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
[926]	1442	for i in range(5):
	1443	dMdpk[i] += 100.cw[iBeg:iFin]refl[13]refl[9]dMdipk[i]
	1444	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
	1445	if Ka2:
	1446	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
	1447	iBeg2 = np.searchsorted(x,pos2-fmin)
	1448	iFin2 = np.searchsorted(x,pos2+fmax)
	1449	if iBeg2-iFin2:
	1450	lenBF2 = iFin2-iBeg2
	1451	dMdpk2 = np.zeros(shape=(6,lenBF2))
[1025]	1452	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
[926]	1453	for i in range(5):
	1454	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[13]refl[9]kRatio*dMdipk2[i]
	1455	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[13]*dMdipk2[0]
	1456	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
	1457	if Phase['General'].get('doPawley'):
	1458	dMdpw = np.zeros(len(x))
	1459	try:
	1460	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
	1461	idx = varylist.index(pIdx)
	1462	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
	1463	if Ka2:
	1464	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
	1465	dMdv[idx] = dMdpw
	1466	except: # ValueError:
	1467	pass
	1468	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
	1469	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
	1470	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
	1471	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
	1472	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
	1473	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
	1474	hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
	1475	for name in names:
	1476	item = names[name]
	1477	if name in varylist:
	1478	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
	1479	if Ka2:
	1480	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
	1481	elif name in dependentVars:
	1482	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
	1483	if Ka2:
	1484	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
	1485	for iPO in dIdPO:
	1486	if iPO in varylist:
	1487	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
	1488	if Ka2:
	1489	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
	1490	elif iPO in dependentVars:
	1491	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
	1492	if Ka2:
	1493	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
	1494	for i,name in enumerate(['omega','chi','phi']):
	1495	aname = pfx+'SH '+name
	1496	if aname in varylist:
	1497	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
	1498	if Ka2:
	1499	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
	1500	elif aname in dependentVars:
	1501	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
	1502	if Ka2:
	1503	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
	1504	for iSH in dFdODF:
	1505	if iSH in varylist:
	1506	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
	1507	if Ka2:
	1508	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
	1509	elif iSH in dependentVars:
	1510	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
	1511	if Ka2:
	1512	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
	1513	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
	1514	for name,dpdA in cellDervNames:
	1515	if name in varylist:
	1516	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
	1517	if Ka2:
	1518	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
	1519	elif name in dependentVars:
	1520	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
	1521	if Ka2:
	1522	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
	1523	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
	1524	for name in dDijDict:
	1525	if name in varylist:
	1526	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
	1527	if Ka2:
	1528	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
	1529	elif name in dependentVars:
	1530	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
	1531	if Ka2:
	1532	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
	1533	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
	1534	for name in gamDict:
	1535	if name in varylist:
	1536	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
	1537	if Ka2:
	1538	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
	1539	elif name in dependentVars:
	1540	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
	1541	if Ka2:
	1542	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
	1543	for name in sigDict:
	1544	if name in varylist:
	1545	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
	1546	if Ka2:
	1547	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
	1548	elif name in dependentVars:
	1549	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
	1550	if Ka2:
	1551	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
	1552	for name in ['BabA','BabU']:
	1553	if phfx+name in varylist:
	1554	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]dervDict['int']cw[iBeg:iFin]
	1555	if Ka2:
	1556	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]dervDict2['int']cw[iBeg2:iFin2]
	1557	elif phfx+name in dependentVars:
	1558	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]dervDict['int']cw[iBeg:iFin]
	1559	if Ka2:
	1560	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]dervDict2['int']cw[iBeg2:iFin2]
	1561	elif 'T' in calcControls[hfx+'histType']:
	1562	print 'TOF Undefined at present'
	1563	raise Exception #no TOF yet
[934]	1564	if not Phase['General'].get('doPawley'):
	1565	#do atom derivatives - for RB,F,X & U so far
	1566	corr = dervDict['int']/refl[9]
	1567	if Ka2:
	1568	corr2 = dervDict2['int']/refl[9]
	1569	for name in varylist+dependentVars:
	1570	if '::RBV;' in name:
	1571	pass
	1572	else:
	1573	try:
	1574	aname = name.split(pfx)[1][:2]
	1575	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
	1576	except IndexError:
	1577	continue
	1578	if name in varylist:
	1579	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
	1580	if Ka2:
	1581	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
	1582	elif name in dependentVars:
	1583	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
	1584	if Ka2:
	1585	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
[926]	1586	# print 'profile derv time: %.3fs'%(time.time()-time0)
	1587	# now process derivatives in constraints
	1588	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
	1589	return dMdv
	1590
	1591	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
[990]	1592	'''Loop over histograms and compute derivatives of the fitting
	1593	model (M) with respect to all parameters. Results are returned in
	1594	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
	1595	n is the number of parameters and m is the number of data
	1596	points. This can exceed memory when m gets large. This routine is
	1597	used when refinement derivatives are selected as "analtytic
	1598	Jacobian" in Controls.
	1599
	1600	:returns: Jacobian numpy.array dMdv for all histograms concatinated
	1601	'''
[926]	1602	parmDict.update(zip(varylist,values))
	1603	G2mv.Dict2Map(parmDict,varylist)
	1604	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	1605	nvar = len(varylist)
	1606	dMdv = np.empty(0)
	1607	histoList = Histograms.keys()
	1608	histoList.sort()
	1609	for histogram in histoList:
	1610	if 'PWDR' in histogram[:4]:
	1611	Histogram = Histograms[histogram]
	1612	hId = Histogram['hId']
	1613	hfx = ':%d:'%(hId)
	1614	wtFactor = calcControls[hfx+'wtFactor']
	1615	Limits = calcControls[hfx+'Limits']
	1616	x,y,w,yc,yb,yd = Histogram['Data']
	1617	W = wtFactor*w
	1618	xB = np.searchsorted(x,Limits[0])
	1619	xF = np.searchsorted(x,Limits[1])
	1620	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
	1621	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
	1622	elif 'HKLF' in histogram[:4]:
	1623	Histogram = Histograms[histogram]
[1030]	1624	nobs = Histogram['Residuals']['Nobs']
[926]	1625	phase = Histogram['Reflection Lists']
	1626	Phase = Phases[phase]
	1627	hId = Histogram['hId']
	1628	hfx = ':%d:'%(hId)
	1629	wtFactor = calcControls[hfx+'wtFactor']
	1630	pfx = '%d::'%(Phase['pId'])
	1631	phfx = '%d:%d:'%(Phase['pId'],hId)
	1632	SGData = Phase['General']['SGData']
	1633	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	1634	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	1635	refList = Histogram['Data']
	1636	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
	1637	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
	1638	dMdvh = np.zeros((len(varylist),len(refList)))
[989]	1639	dependentVars = G2mv.GetDependentVars()
	1640	depDerivDict = {}
	1641	for j in dependentVars:
	1642	depDerivDict[j] = np.zeros(shape=(len(refList)))
[926]	1643	if calcControls['F**2']:
	1644	for iref,ref in enumerate(refList):
	1645	if ref[6] > 0:
	1646	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
	1647	w = 1.0/ref[6]
	1648	if w*ref[5] >= calcControls['minF/sig']:
	1649	for j,var in enumerate(varylist):
	1650	if var in dFdvDict:
	1651	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*dervCor
[989]	1652	for var in dependentVars:
	1653	if var in dFdvDict:
	1654	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*dervCor
[926]	1655	if phfx+'Scale' in varylist:
	1656	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
[989]	1657	elif phfx+'Scale' in dependentVars:
	1658	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
[926]	1659	for item in ['Ep','Es','Eg']:
	1660	if phfx+item in varylist:
	1661	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[989]	1662	elif phfx+item in dependentVars:
	1663	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[926]	1664	for item in ['BabA','BabU']:
	1665	if phfx+item in varylist:
	1666	dMdvh[varylist.index(phfx+item)][iref] = wdervCordFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
[989]	1667	elif phfx+item in dependentVars:
	1668	depDerivDict[phfx+item][iref] = wdervCordFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
[926]	1669	else:
	1670	for iref,ref in enumerate(refList):
	1671	if ref[5] > 0.:
	1672	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
	1673	Fo = np.sqrt(ref[5])
	1674	Fc = np.sqrt(ref[7])
	1675	w = 1.0/ref[6]
	1676	if 2.0Fow*Fo >= calcControls['minF/sig']:
	1677	for j,var in enumerate(varylist):
	1678	if var in dFdvDict:
	1679	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
[989]	1680	for var in dependentVars:
	1681	if var in dFdvDict:
	1682	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
[926]	1683	if phfx+'Scale' in varylist:
[989]	1684	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
	1685	elif phfx+'Scale' in dependentVars:
	1686	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
[926]	1687	for item in ['Ep','Es','Eg']:
	1688	if phfx+item in varylist:
	1689	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[989]	1690	elif phfx+item in dependentVars:
	1691	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[926]	1692	for item in ['BabA','BabU']:
	1693	if phfx+item in varylist:
	1694	dMdvh[varylist.index(phfx+item)][iref] = wdervCordFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
[989]	1695	elif phfx+item in dependentVars:
	1696	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
	1697	# now process derivatives in constraints
	1698	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
[926]	1699	else:
	1700	continue #skip non-histogram entries
	1701	if len(dMdv):
	1702	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
	1703	else:
	1704	dMdv = np.sqrt(wtFactor)*dMdvh
	1705
	1706	pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
	1707	if np.any(pVals):
	1708	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
	1709	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
	1710
	1711	return dMdv
	1712
	1713	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
[990]	1714	'''Loop over histograms and compute derivatives of the fitting
	1715	model (M) with respect to all parameters. For each histogram, the
	1716	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
	1717	number of parameters and m is the number of data points *in the
	1718	histogram*. The (n by n) Hessian is computed from each Jacobian
	1719	and it is returned. This routine is used when refinement
	1720	derivatives are selected as "analtytic Hessian" in Controls.
	1721
	1722	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
	1723	'''
[939]	1724	'Needs a doc string'
[926]	1725	parmDict.update(zip(varylist,values))
	1726	G2mv.Dict2Map(parmDict,varylist)
	1727	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	1728	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
	1729	nvar = len(varylist)
	1730	Hess = np.empty(0)
	1731	histoList = Histograms.keys()
	1732	histoList.sort()
	1733	for histogram in histoList:
	1734	if 'PWDR' in histogram[:4]:
	1735	Histogram = Histograms[histogram]
	1736	hId = Histogram['hId']
	1737	hfx = ':%d:'%(hId)
	1738	wtFactor = calcControls[hfx+'wtFactor']
	1739	Limits = calcControls[hfx+'Limits']
	1740	x,y,w,yc,yb,yd = Histogram['Data']
	1741	W = wtFactor*w
	1742	dy = y-yc
	1743	xB = np.searchsorted(x,Limits[0])
	1744	xF = np.searchsorted(x,Limits[1])
	1745	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
	1746	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
[1017]	1747	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
[926]	1748	Dy = dy[xB:xF][np.newaxis,:]
	1749	dMdvh *= Wt
	1750	if dlg:
[1030]	1751	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
[926]	1752	if len(Hess):
	1753	Hess += np.inner(dMdvh,dMdvh)
	1754	dMdvh = WtDy
	1755	Vec += np.sum(dMdvh,axis=1)
	1756	else:
	1757	Hess = np.inner(dMdvh,dMdvh)
	1758	dMdvh = WtDy
	1759	Vec = np.sum(dMdvh,axis=1)
	1760	elif 'HKLF' in histogram[:4]:
	1761	Histogram = Histograms[histogram]
[1030]	1762	nobs = Histogram['Residuals']['Nobs']
[926]	1763	phase = Histogram['Reflection Lists']
	1764	Phase = Phases[phase]
	1765	hId = Histogram['hId']
	1766	hfx = ':%d:'%(hId)
	1767	wtFactor = calcControls[hfx+'wtFactor']
	1768	pfx = '%d::'%(Phase['pId'])
	1769	phfx = '%d:%d:'%(Phase['pId'],hId)
	1770	SGData = Phase['General']['SGData']
	1771	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	1772	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	1773	refList = Histogram['Data']
	1774	time0 = time.time()
	1775	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict) #accurate for powders!
	1776	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
	1777	dMdvh = np.zeros((len(varylist),len(refList)))
[989]	1778	dependentVars = G2mv.GetDependentVars()
	1779	depDerivDict = {}
	1780	for j in dependentVars:
	1781	depDerivDict[j] = np.zeros(shape=(len(refList)))
[926]	1782	wdf = np.zeros(len(refList))
	1783	if calcControls['F**2']:
	1784	for iref,ref in enumerate(refList):
	1785	if ref[6] > 0:
	1786	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
	1787	w = 1.0/ref[6]
	1788	if w*ref[5] >= calcControls['minF/sig']:
	1789	wdf[iref] = w*(ref[5]-ref[7])
	1790	for j,var in enumerate(varylist):
	1791	if var in dFdvDict:
	1792	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
[989]	1793	for var in dependentVars:
	1794	if var in dFdvDict:
	1795	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
[926]	1796	if phfx+'Scale' in varylist:
	1797	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
[989]	1798	elif phfx+'Scale' in dependentVars:
	1799	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
[926]	1800	for item in ['Ep','Es','Eg']:
	1801	if phfx+item in varylist:
	1802	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[989]	1803	elif phfx+item in dependentVars:
	1804	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[926]	1805	for item in ['BabA','BabU']:
	1806	if phfx+item in varylist:
	1807	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
[989]	1808	elif phfx+item in dependentVars:
	1809	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
[926]	1810	else:
	1811	for iref,ref in enumerate(refList):
	1812	if ref[5] > 0.:
	1813	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
	1814	Fo = np.sqrt(ref[5])
	1815	Fc = np.sqrt(ref[7])
	1816	w = 1.0/ref[6]
	1817	if 2.0Fow*Fo >= calcControls['minF/sig']:
	1818	wdf[iref] = 2.0Fow*(Fo-Fc)
	1819	for j,var in enumerate(varylist):
	1820	if var in dFdvDict:
	1821	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
[989]	1822	for var in dependentVars:
	1823	if var in dFdvDict:
	1824	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
[926]	1825	if phfx+'Scale' in varylist:
[989]	1826	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
	1827	elif phfx+'Scale' in dependentVars:
	1828	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
[926]	1829	for item in ['Ep','Es','Eg']:
	1830	if phfx+item in varylist:
	1831	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[989]	1832	elif phfx+item in dependentVars:
	1833	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
[926]	1834	for item in ['BabA','BabU']:
	1835	if phfx+item in varylist:
	1836	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
[989]	1837	elif phfx+item in dependentVars:
	1838	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
	1839	# now process derivatives in constraints
	1840	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
	1841
[926]	1842	if dlg:
[1030]	1843	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
[926]	1844	if len(Hess):
	1845	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
	1846	Hess += wtFactor*np.inner(dMdvh,dMdvh)
	1847	else:
	1848	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
	1849	Hess = wtFactor*np.inner(dMdvh,dMdvh)
	1850	else:
	1851	continue #skip non-histogram entries
	1852	pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
	1853	if np.any(pVals):
	1854	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
	1855	Vec += np.sum(dpdvpWtpVals,axis=1)
	1856	Hess += np.inner(dpdv*pWt,dpdv)
	1857	return Vec,Hess
	1858
	1859	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
[939]	1860	'Needs a doc string'
[926]	1861	parmDict.update(zip(varylist,values))
	1862	Values2Dict(parmDict, varylist, values)
	1863	G2mv.Dict2Map(parmDict,varylist)
	1864	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
	1865	M = np.empty(0)
	1866	SumwYo = 0
	1867	Nobs = 0
	1868	ApplyRBModels(parmDict,Phases,rigidbodyDict)
	1869	histoList = Histograms.keys()
	1870	histoList.sort()
	1871	for histogram in histoList:
	1872	if 'PWDR' in histogram[:4]:
	1873	Histogram = Histograms[histogram]
	1874	hId = Histogram['hId']
	1875	hfx = ':%d:'%(hId)
	1876	wtFactor = calcControls[hfx+'wtFactor']
	1877	Limits = calcControls[hfx+'Limits']
	1878	x,y,w,yc,yb,yd = Histogram['Data']
	1879	yc *= 0.0 #zero full calcd profiles
	1880	yb *= 0.0
	1881	yd *= 0.0
	1882	xB = np.searchsorted(x,Limits[0])
	1883	xF = np.searchsorted(x,Limits[1])
	1884	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
	1885	varylist,Histogram,Phases,calcControls,pawleyLookup)
	1886	yc[xB:xF] += yb[xB:xF]
[1087]	1887	if not np.any(y): #fill dummy data
	1888	rv = st.poisson(yc[xB:xF])
	1889	y[xB:xF] = rv.rvs()
[1088]	1890	Z = np.ones_like(yc[xB:xF])
	1891	Z[1::2] *= -1
	1892	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
[1087]	1893	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
[926]	1894	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
[1087]	1895	W = wtFactor*w
[1017]	1896	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
[1087]	1897	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
	1898	Nobs += Histogram['Residuals']['Nobs']
	1899	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
	1900	SumwYo += Histogram['Residuals']['sumwYo']
[1038]	1901	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
[1030]	1902	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
[1038]	1903	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
	1904	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
	1905	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.))
	1906	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.))
	1907	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
	1908	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
	1909	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
[926]	1910	if dlg:
[1030]	1911	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
[926]	1912	M = np.concatenate((M,wdy))
	1913	#end of PWDR processing
	1914	elif 'HKLF' in histogram[:4]:
	1915	Histogram = Histograms[histogram]
[1043]	1916	Histogram['Residuals'] = {}
[926]	1917	phase = Histogram['Reflection Lists']
	1918	Phase = Phases[phase]
	1919	hId = Histogram['hId']
	1920	hfx = ':%d:'%(hId)
	1921	wtFactor = calcControls[hfx+'wtFactor']
	1922	pfx = '%d::'%(Phase['pId'])
	1923	phfx = '%d:%d:'%(Phase['pId'],hId)
	1924	SGData = Phase['General']['SGData']
	1925	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
	1926	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
	1927	refList = Histogram['Data']
	1928	time0 = time.time()
	1929	StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
	1930	# print 'sf-calc time: %.3f'%(time.time()-time0)
	1931	df = np.zeros(len(refList))
	1932	sumwYo = 0
	1933	sumFo = 0
	1934	sumFo2 = 0
	1935	sumdF = 0
	1936	sumdF2 = 0
	1937	nobs = 0
	1938	if calcControls['F**2']:
	1939	for i,ref in enumerate(refList):
	1940	if ref[6] > 0:
	1941	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
	1942	w = 1.0/ref[6]
	1943	ref[7] = parmDict[phfx+'Scale']*ref[9]
	1944	ref[7] *= ref[13] #correct Fc^2 for extinction
	1945	ref[8] = ref[5]/parmDict[phfx+'Scale']
	1946	if w*ref[5] >= calcControls['minF/sig']:
	1947	sumFo2 += ref[5]
	1948	Fo = np.sqrt(ref[5])
	1949	sumFo += Fo
	1950	sumFo2 += ref[5]
	1951	sumdF += abs(Fo-np.sqrt(ref[7]))
	1952	sumdF2 += abs(ref[5]-ref[7])
	1953	nobs += 1
	1954	df[i] = -w*(ref[5]-ref[7])
	1955	sumwYo += (wref[5])*2
	1956	else:
	1957	for i,ref in enumerate(refList):
	1958	if ref[5] > 0.:
	1959	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
	1960	ref[7] = parmDict[phfx+'Scale']*ref[9]
	1961	ref[7] *= ref[13] #correct Fc^2 for extinction
	1962	Fo = np.sqrt(ref[5])
	1963	Fc = np.sqrt(ref[7])
	1964	w = 2.0*Fo/ref[6]
	1965	if w*Fo >= calcControls['minF/sig']:
	1966	sumFo += Fo
	1967	sumFo2 += ref[5]
	1968	sumdF += abs(Fo-Fc)
	1969	sumdF2 += abs(ref[5]-ref[7])
	1970	nobs += 1
	1971	df[i] = -w*(Fo-Fc)
	1972	sumwYo += (wFo)*2
[1030]	1973	Histogram['Residuals']['Nobs'] = nobs
	1974	Histogram['Residuals']['sumwYo'] = sumwYo
[926]	1975	SumwYo += sumwYo
[1030]	1976	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['Residuals']['sumwYo'])100.)
	1977	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
	1978	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
	1979	Histogram['Residuals'][phfx+'Nref'] = nobs
[926]	1980	Nobs += nobs
	1981	if dlg:
[1030]	1982	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
[926]	1983	M = np.concatenate((M,wtFactor*df))
	1984	# end of HKLF processing
	1985	Histograms['sumwYo'] = SumwYo
	1986	Histograms['Nobs'] = Nobs
	1987	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
	1988	if dlg:
	1989	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
	1990	if not GoOn:
	1991	parmDict['saved values'] = values
	1992	dlg.Destroy()
	1993	raise Exception #Abort!!
	1994	pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
	1995	if np.any(pVals):
	1996	pSum = np.sum(pWtpVals*2)
	1997	print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
	1998	Nobs += len(pVals)
	1999	M = np.concatenate((M,np.sqrt(pWt)*pVals))
	2000	return M
	2001

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

Download in other formats: