Changeset 1950
- Timestamp:
- Aug 4, 2015 3:51:43 PM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIspc.py
r1948 r1950 1430 1430 return CSuinel[indx[1]] 1431 1431 1432 #def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False): 1433 # 1434 # def fracCrenel(tau,Toff,Twid): 1435 # Tau = (tau-Toff)%1. 1436 # A = np.where(Tau<Twid,1.,0.) 1437 # return A 1438 # 1439 # def fracFourier(tau,nH,fsin,fcos): 1440 # SA = np.sin(2.*nH*np.pi*tau) 1441 # CB = np.cos(2.*nH*np.pi*tau) 1442 # A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis] 1443 # B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis] 1444 # return A+B 1445 # 1446 # def posFourier(tau,nH,psin,pcos): 1447 # SA = np.sin(2*nH*np.pi*tau) 1448 # CB = np.cos(2*nH*np.pi*tau) 1449 # A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis] 1450 # B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis] 1451 # return A+B 1452 # 1453 # def posSawtooth(tau,Toff,slopes): 1454 # Tau = (tau-Toff[:,np.newaxis])%1. 1455 # A = slopes[:,:,np.newaxis]*Tau 1456 # return A 1457 # 1458 # def posZigZag(tau,Toff,slopes): 1459 # Tau = (tau-Toff[:,np.newaxis])%1. 1460 # A = np.where(Tau <= 0.5,slopes[:,:,np.newaxis]*Tau,slopes[:,:,np.newaxis]*(1.-Tau)) 1461 # return A 1462 # 1463 # print 'super space group: ',SSGData['SSpGrp'] 1464 # CSI = {'Sfrac':[[[1,0],[2,0]],[[1.,0.],[1.,0.]]], 1465 # 'Spos':[[[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]], 1466 # [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]], #sin & cos 1467 # 'Sadp':[[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 1468 # [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]], 1469 # [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], 1470 # [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]], 1471 # 'Smag':[[[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]], 1472 # [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],} 1473 # xyz = np.array(XYZ)%1. 1474 # xyzt = np.array(XYZ+[0,])%1. 1475 # SGOps = copy.deepcopy(SGData['SGOps']) 1476 # siteSym = SytSym(XYZ,SGData)[0].strip() 1477 # print 'siteSym: ',siteSym 1478 # if siteSym == '1': #"1" site symmetry 1479 # if debug: 1480 # return CSI,None,None,None,None 1481 # else: 1482 # return CSI 1483 # elif siteSym == '-1': #"-1" site symmetry 1484 # CSI['Sfrac'][0] = [[1,0],[0,0]] 1485 # CSI['Spos'][0] = [[1,0,0],[2,0,0],[3,0,0], [0,0,0],[0,0,0],[0,0,0]] 1486 # CSI['Sadp'][0] = [[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0], 1487 # [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]] 1488 # if debug: 1489 # return CSI,None,None,None,None 1490 # else: 1491 # return CSI 1492 # SSGOps = copy.deepcopy(SSGData['SSGOps']) 1493 # #expand ops to include inversions if any 1494 # if SGData['SGInv']: 1495 # for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']): 1496 # SGOps.append([-op[0],-op[1]%1.]) 1497 # SSGOps.append([-sop[0],-sop[1]%1.]) 1498 # #build set of sym ops around special poasition 1499 # SSop = [] 1500 # Sop = [] 1501 # Sdtau = [] 1502 # for iop,Op in enumerate(SGOps): 1503 # nxyz = (np.inner(Op[0],xyz)+Op[1])%1. 1504 # if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z': 1505 # SSop.append(SSGOps[iop]) 1506 # Sop.append(SGOps[iop]) 1507 # ssopinv = nl.inv(SSGOps[iop][0]) 1508 # mst = ssopinv[3][:3] 1509 # epsinv = ssopinv[3][3] 1510 # Sdtau.append(np.sum(mst*(XYZ-SGOps[iop][1])-epsinv*SSGOps[iop][1][3])) 1511 # SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order 1512 # OpText = [MT2text(s).replace(' ','') for s in Sop] #debug? 1513 # SSOpText = [SSMT2text(ss).replace(' ','') for ss in SSop] #debug? 1514 # print 'special pos super operators: ',SSOpText 1515 # #setup displacement arrays 1516 # tau = np.linspace(0,1,49,True) 1517 # delt2 = np.eye(2)*0.001 1518 # delt4 = np.eye(4)*0.001 1519 # delt6 = np.eye(6)*0.001 1520 # delt12 = np.eye(12)*0.0001 1521 # #make modulation arrays - one parameter at a time 1522 # #site fractions 1523 # CSI['Sfrac'] = [np.zeros((2),dtype='i'),np.ones(2)] 1524 # if 'Crenel' in waveType: 1525 # dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze() 1526 # else: 1527 # dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze() 1528 # dFT = np.zeros_like(dF) 1529 # #positions 1530 # if 'Fourier' in waveType: 1531 # dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis] 1532 # #3x6x12 modulated position array (X,Spos,tau)& force positive 1533 # CSI['Spos'] = [np.zeros((6,3),dtype='i'),np.zeros((6,3))] 1534 # elif waveType == 'Sawtooth': 1535 # CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1536 # elif waveType == 'ZigZag': 1537 # CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1538 # #anisotropic thermal motion 1539 # dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array 1540 # CSI['Sadp'] = [np.zeros((12,3),dtype='i'),np.zeros((12,3))] 1541 # 1542 # FSC = np.ones(2,dtype='i') 1543 # VFSC = np.ones(2) 1544 # XSC = np.ones(6,dtype='i') 1545 # USC = np.ones(12,dtype='i') 1546 # dFTP = [] 1547 # dXTP = [] 1548 # dUTP = [] 1549 # for i in SdIndx: 1550 # sop = Sop[i] 1551 # ssop = SSop[i] 1552 # fsc = np.ones(2,dtype='i') 1553 # xsc = np.ones(6,dtype='i') 1554 # ssopinv = nl.inv(ssop[0]) 1555 # mst = ssopinv[3][:3] 1556 # epsinv = ssopinv[3][3] 1557 # sdet = nl.det(sop[0]) 1558 # ssdet = nl.det(ssop[0]) 1559 # dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3] 1560 # dT = 1.0 1561 # if np.any(dtau%.5): 1562 # dT = np.tan(np.pi*np.sum(dtau%.5)) 1563 # tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]) 1564 # if waveType == 'Fourier': 1565 # dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis] 1566 # elif waveType == 'Sawtooth': 1567 # dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1568 # elif waveType == 'ZigZag': 1569 # dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1570 # dXT = np.inner(sop[0],dXT.T) 1571 # dXT = np.swapaxes(dXT,1,2) 1572 # dXT[:,:3,:] *= ssdet 1573 # dXTP.append(dXT) 1574 # if waveType == 'Fourier': 1575 # if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1576 # CSI['Spos'] = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]], 1577 # [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]] 1578 # if '(x)' in siteSym: 1579 # CSI['Spos'][1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.] 1580 # if 'm' in siteSym and len(SdIndx) == 1: 1581 # CSI['Spos'][1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.] 1582 # elif '(y)' in siteSym: 1583 # CSI['Spos'][1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1584 # if 'm' in siteSym and len(SdIndx) == 1: 1585 # CSI['Spos'][1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1586 # elif '(z)' in siteSym: 1587 # CSI['Spos'][1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1588 # if 'm' in siteSym and len(SdIndx) == 1: 1589 # CSI['Spos'][1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1590 # for i in range(3): 1591 # if not XSC[i]: 1592 # CSI['Spos'][0][i] = [0,0,0] 1593 # CSI['Spos'][1][i] = [0.,0.,0.] 1594 # CSI['Spos'][0][i+3] = [0,0,0] 1595 # CSI['Spos'][1][i+3] = [0.,0.,0.] 1596 # else: 1597 # for i in range(3): 1598 # if np.allclose(dX[i,i,:],dXT[i,i,:]*sdet): 1599 # xsc[i] = 1 1600 # else: 1601 # xsc[i] = 0 1602 # if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]): 1603 # xsc[i+3] = 1 1604 # else: 1605 # xsc[i+3] = 0 1606 # XSC &= xsc 1607 # 1608 # fsc = np.ones(2,dtype='i') 1609 # if 'Crenel' in waveType: 1610 # dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze() 1611 # fsc = [1,1] 1612 # else: 1613 # dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze() 1614 # dFT = nl.det(sop[0])*dFT 1615 # dFT = dFT[:,np.argsort(tauT)] 1616 # dFT[0] *= ssdet 1617 # dFT[1] *= sdet 1618 # dFTP.append(dFT) 1619 # 1620 # if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1621 # fsc = [1,1] 1622 # CSI['Sfrac'] = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]] 1623 # for i in range(2): 1624 # if not FSC[i]: 1625 # CSI['Sfrac'][0][i] = [0,0] 1626 # CSI['Sfrac'][1][i] = [0.,0.] 1627 # else: 1628 # for i in range(2): 1629 # if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6): 1630 # fsc[i] = 1 1631 # else: 1632 # fsc[i] = 0 1633 # FSC &= fsc 1634 # 1635 # usc = np.ones(12,dtype='i') 1636 # dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array 1637 # dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3 1638 # dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) 1639 # dUT = np.array(U2Uij(dUijT)) 1640 # dUT = dUT[:,:,np.argsort(tauT)] 1641 # dUTP.append(dUT) 1642 # if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1643 # CSI['Sadp'] = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 1644 # [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], 1645 # [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], 1646 # [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]] 1647 # if '(x)' in siteSym: 1648 # CSI['Sadp'][1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1649 # elif '(y)' in siteSym: 1650 # CSI['Sadp'][1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1651 # elif '(z)' in siteSym: 1652 # CSI['Sadp'][1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.] 1653 # for i in range(6): 1654 # if not USC[i]: 1655 # CSI['Sadp'][0][i] = [0,0,0] 1656 # CSI['Sadp'][1][i] = [0.,0.,0.] 1657 # CSI['Sadp'][0][i+6] = [0,0,0] 1658 # CSI['Sadp'][1][i+6] = [0.,0.,0.] 1659 # else: 1660 # 1661 # for i in range(6): 1662 # if np.allclose(dU[i,i,:],dUT[i,i,:]*sdet): 1663 # usc[i] = 1 1664 # else: 1665 # usc[i] = 0 1666 # if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): 1667 # usc[i+6] = 1 1668 # else: 1669 # usc[i+6] = 0 1670 # if '4/m' in siteSym and np.any(dUT[0,1,:]): 1671 # CSI['Sadp'][0][6:8] = [[12,0,0],[12,0,0]] 1672 # if ssop[1][3]: 1673 # CSI['Sadp'][1][6:8] = [[1.,0.,0.],[-1.,0.,0.]] 1674 # usc[9] = 1 1675 # else: 1676 # CSI['Sadp'][1][6:8] = [[1.,0.,0.],[1.,0.,0.]] 1677 # usc[9] = 0 1678 # elif '4' in siteSym and np.any(dUT[0,1,:]): 1679 # CSI['Sadp'][0][6:8] = [[12,0,0],[12,0,0]] 1680 # CSI['Sadp'][0][:2] = [[11,0,0],[11,0,0]] 1681 # if ssop[1][3]: 1682 # CSI['Sadp'][1][:2] = [[1.,0.,0.],[-1.,0.,0.]] 1683 # CSI['Sadp'][1][6:8] = [[1.,0.,0.],[-1.,0.,0.]] 1684 # usc[3] = 1 1685 # usc[9] = 1 1686 # else: 1687 # CSI['Sadp'][1][:2] = [[1.,0.,0.],[1.,0.,0.]] 1688 # CSI['Sadp'][1][6:8] = [[1.,0.,0.],[1.,0.,0.]] 1689 # usc[3] = 0 1690 # usc[9] = 0 1691 # print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc 1692 # USC &= usc 1693 # print dtau 1694 # print FSC 1695 # print XSC 1696 # print USC 1697 # if not np.any(dtau%.5): 1698 # n = -1 1699 # for i,U in enumerate(USC): 1700 # if U: 1701 # n += 1 1702 # CSI['Sadp'][0][i][0] = n+1 1703 # CSI['Sadp'][1][i][0] = 1.0 1704 # if waveType == 'Fourier': 1705 # n = -1 1706 # for i,X in enumerate(XSC): 1707 # if X: 1708 # n += 1 1709 # CSI['Spos'][0][i][0] = n+1 1710 # CSI['Spos'][1][i][0] = 1.0 1711 # n = -1 1712 # for i,F in enumerate(FSC): 1713 # if F: 1714 # n += 1 1715 # CSI['Sfrac'][0][i] = n+1 1716 # CSI['Sfrac'][1][i] = 1.0 1717 # else: 1718 # CSI['Sfrac'][0][i] = 0 1719 # CSI['Sfrac'][1][i] = 0. 1720 # if debug: 1721 # return CSI,[tau,tauT],[dF,dFTP],[dX,dXTP],[dU,dUTP] 1722 # else: 1723 # return CSI 1724 # 1432 1725 def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False): 1433 1726 … … 1444 1737 1445 1738 def fracCrenel(tau,Toff,Twid): 1446 Tau = (tau-Toff )%1.1447 A = np.where(Tau<Twid ,1.,0.)1739 Tau = (tau-Toff[:,np.newaxis])%1. 1740 A = np.where(Tau<Twid[:,np.newaxis],1.,0.) 1448 1741 return A 1449 1742 … … 1463 1756 1464 1757 def posSawtooth(tau,Toff,slopes): 1465 Tau = (tau-Toff [:,np.newaxis])%1.1466 A = slopes[:, :,np.newaxis]*Tau1758 Tau = (tau-Toff)%1. 1759 A = slopes[:,np.newaxis]*Tau 1467 1760 return A 1468 1761 1469 1762 def posZigZag(tau,Toff,slopes): 1470 Tau = (tau-Toff [:,np.newaxis])%1.1471 A = np.where(Tau <= 0.5,slopes[:, :,np.newaxis]*Tau,slopes[:,:,np.newaxis]*(1.-Tau))1763 Tau = (tau-Toff)%1. 1764 A = np.where(Tau <= 0.5,slopes[:,np.newaxis]*Tau,slopes[:,np.newaxis]*(1.-Tau)) 1472 1765 return A 1473 1766 … … 1531 1824 1532 1825 def DoXYZ(): 1826 delt4 = np.ones(4)*0.001 1533 1827 delt6 = np.eye(6)*0.001 1534 1828 if 'Fourier' in waveType: … … 1537 1831 CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))] 1538 1832 elif waveType == 'Sawtooth': 1539 CSI = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1833 dX = posSawtooth(tau,delt4[0],delt4[1:]) 1834 CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]), 1835 np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])] 1540 1836 elif waveType == 'ZigZag': 1541 CSI = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1837 dX = posZigZag(tau,delt4[0],delt4[1:]) 1838 CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]), 1839 np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])] 1542 1840 XSC = np.ones(6,dtype='i') 1543 1841 dXTP = [] … … 1561 1859 dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1562 1860 elif waveType == 'ZigZag': 1563 dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1564 dXT = np.inner(sop[0],dXT.T) 1565 dXT = np.swapaxes(dXT,1,2) 1566 dXT[:,:3,:] *= ssdet1861 dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1862 dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6) 1863 dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49) 1864 dXT[:,:3,:] *= (ssdet*sdet) # modify the sin component 1567 1865 dXTP.append(dXT) 1866 # print np.sum(dX**2,axis=-1)*1e5 1867 # print np.sum(dXT**2,axis=-1)*1e5 1568 1868 if waveType == 'Fourier': 1569 1869 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): … … 1586 1886 CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]] 1587 1887 CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]] 1588 elif np.allclose(dX[0,0,:],dXT[0,1,:]*sdet): 1589 if 'xy' in siteSym or '+-0' in siteSym: 1590 CSI[0][3:5] = [[12,0,0],[12,0,0]] 1591 CSI[1][3:5] = [[1.,0,0],[-sdet,0,0]] 1888 elif '(xy)' in siteSym or '(+-0)' in siteSym: 1889 mul = 1 1890 if '(xy)' in siteSym: 1891 mul = -1 1892 if np.allclose(dX[0,0,:],dXT[1,0,:]*mul): 1893 CSI[0][3:5] = [[11,0,0],[11,0,0]] 1894 CSI[1][3:5] = [[1.,0,0],[mul,0,0]] 1592 1895 xsc[3:5] = 0 1896 if np.allclose(dX[0,3,:],dXT[0,4,:]*mul): 1897 CSI[0][:2] = [[12,0,0],[12,0,0]] 1898 CSI[1][:2] = [[1.,0,0],[mul,0,0]] 1899 xsc[:2] = 0 1593 1900 else: 1594 1901 for i in range(3): 1595 if np.allclose(dX[i,i,:],dXT[i,i,:]*sdet): 1596 xsc[i] = 1 1597 else: 1902 if not np.allclose(dX[i,i,:],dXT[i,i,:]): 1598 1903 xsc[i] = 0 1599 if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]): 1600 xsc[i+3] = 1 1601 else: 1904 if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]): 1602 1905 xsc[i+3] = 0 1603 1906 XSC &= xsc 1604 1907 print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc 1605 n = -1 1606 print XSC 1607 for i,X in enumerate(XSC): 1608 if X: 1609 n += 1 1610 CSI[0][i][0] = n+1 1611 CSI[1][i][0] = 1.0 1908 if waveType == 'Fourier': 1909 n = -1 1910 print XSC 1911 for i,X in enumerate(XSC): 1912 if X: 1913 n += 1 1914 CSI[0][i][0] = n+1 1915 CSI[1][i][0] = 1.0 1612 1916 1613 1917 return CSI,dX,dXTP 1614 1918 1615 1919 def DoUij(): 1920 tau = np.linspace(0,1,49,True) 1616 1921 delt12 = np.eye(12)*0.0001 1617 1922 dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array … … 1635 1940 dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array 1636 1941 dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3 1637 dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) 1638 dUT = np.array(U2Uij(dUijT)) 1942 dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49 1943 dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49 1639 1944 dUT = dUT[:,:,np.argsort(tauT)] 1640 1945 dUTP.append(dUT) … … 1666 1971 else: 1667 1972 for i in range(6): 1668 if np.allclose(dU[i,i,:],-dUT[i,i,:]): 1669 usc[i] = 1 1670 else: 1973 if not np.allclose(dU[i,i,:],dUT[i,i,:]*sdet): #sin part 1671 1974 usc[i] = 0 1672 if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): 1673 usc[i+6] = 1 1674 else: 1975 if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part 1675 1976 usc[i+6] = 0 1676 if np.any(dUT[ 0,1,:]):1977 if np.any(dUT[1,0,:]): 1677 1978 if '4/m' in siteSym: 1678 1979 CSI[0][6:8] = [[12,0,0],[12,0,0]] … … 1702 2003 elif 'xy' in siteSym or '+-0' in siteSym: 1703 2004 if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet): 1704 print np.allclose(dU[0,0,:],dUT[0,1,:]*sdet),sdet1705 2005 CSI[0][4:6] = [[12,0,0],[12,0,0]] 1706 2006 CSI[0][6:8] = [[11,0,0],[11,0,0]] … … 1712 2012 print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc 1713 2013 USC &= usc 2014 print USC 1714 2015 if not np.any(dtau%.5): 1715 2016 n = -1 … … 1790 2091 CSI['Sadp'][0] = orderParms(CSI['Sadp'][0]) 1791 2092 if debug: 1792 return CSI, [tau,tau],[dF,dFTP],[dX,dXTP],[dU,dUTP]2093 return CSI,tau,[dF,dFTP],[dX,dXTP],[dU,dUTP] 1793 2094 else: 1794 2095 return CSI
Note: See TracChangeset
for help on using the changeset viewer.