Changeset 1950


Ignore:
Timestamp:
Aug 4, 2015 3:51:43 PM (6 years ago)
Author:
vondreele
Message:

more work on SS special pos.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIspc.py

    r1948 r1950  
    14301430    return CSuinel[indx[1]]
    14311431   
     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#
    14321725def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False):
    14331726   
     
    14441737   
    14451738    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.)
    14481741        return A
    14491742       
     
    14631756
    14641757    def posSawtooth(tau,Toff,slopes):
    1465         Tau = (tau-Toff[:,np.newaxis])%1.
    1466         A = slopes[:,:,np.newaxis]*Tau
     1758        Tau = (tau-Toff)%1.
     1759        A = slopes[:,np.newaxis]*Tau
    14671760        return A
    14681761   
    14691762    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))
    14721765        return A
    14731766       
     
    15311824       
    15321825    def DoXYZ():
     1826        delt4 = np.ones(4)*0.001
    15331827        delt6 = np.eye(6)*0.001
    15341828        if 'Fourier' in waveType:
     
    15371831            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
    15381832        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]])]
    15401836        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]])]
    15421840        XSC = np.ones(6,dtype='i')
    15431841        dXTP = []
     
    15611859                dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
    15621860            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,:] *= ssdet
     1861                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
    15671865            dXTP.append(dXT)
     1866#            print np.sum(dX**2,axis=-1)*1e5
     1867#            print np.sum(dXT**2,axis=-1)*1e5
    15681868            if waveType == 'Fourier':
    15691869                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
     
    15861886                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
    15871887                        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]]
    15921895                        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
    15931900                else:
    15941901                    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,:]):
    15981903                            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,:]):
    16021905                            xsc[i+3] = 0
    16031906                XSC &= xsc
    16041907                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
    16121916       
    16131917        return CSI,dX,dXTP
    16141918       
    16151919    def DoUij():
     1920        tau = np.linspace(0,1,49,True)
    16161921        delt12 = np.eye(12)*0.0001
    16171922        dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
     
    16351940            dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
    16361941            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
    16391944            dUT = dUT[:,:,np.argsort(tauT)]
    16401945            dUTP.append(dUT)
     
    16661971            else:                       
    16671972                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
    16711974                        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
    16751976                        usc[i+6] = 0
    1676                 if np.any(dUT[0,1,:]):
     1977                if np.any(dUT[1,0,:]):
    16771978                    if '4/m' in siteSym:
    16781979                        CSI[0][6:8] = [[12,0,0],[12,0,0]]
     
    17022003                    elif 'xy' in siteSym or '+-0' in siteSym:
    17032004                        if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
    1704                             print np.allclose(dU[0,0,:],dUT[0,1,:]*sdet),sdet
    17052005                            CSI[0][4:6] = [[12,0,0],[12,0,0]]
    17062006                            CSI[0][6:8] = [[11,0,0],[11,0,0]]
     
    17122012                print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc
    17132013            USC &= usc
     2014        print USC
    17142015        if not np.any(dtau%.5):
    17152016            n = -1
     
    17902091    CSI['Sadp'][0] = orderParms(CSI['Sadp'][0])           
    17912092    if debug:
    1792         return CSI,[tau,tau],[dF,dFTP],[dX,dXTP],[dU,dUTP]
     2093        return CSI,tau,[dF,dFTP],[dX,dXTP],[dU,dUTP]
    17932094    else:
    17942095        return CSI
Note: See TracChangeset for help on using the changeset viewer.