Changeset 379


Ignore:
Timestamp:
Sep 22, 2011 10:11:18 AM (11 years ago)
Author:
vondreele
Message:

finish spherical harmonics preferred orientation

Location:
trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIlattice.py

    r378 r379  
    788788   
    789789def SamAng(Tth,Gangls,Sangl,IFCoup):
     790    """Compute sample orientation angles vs laboratory coord. system
     791    input:
     792        Tth:        Signed theta                                   
     793        Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
     794        Sangl:      Sample angle zeros om-0, chi-0, phi-0         
     795        IFCoup:     =.TRUE. if omega & 2-theta coupled in CW scan
     796    returns: 
     797        psi,gam:    Sample odf angles                             
     798        dPSdA,dGMdA:    Angle zero derivatives
     799    """                         
     800   
     801    rpd = math.pi/180.
    790802    if IFCoup:
    791803        GSomeg = sind(Gangls[2]+Tth)
    … …  
    817829    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
    818830    psi = acosd(BC)
     831   
     832    BD = 1.0-BC**2
     833    if BD > 0.:
     834        C = rpd/math.sqrt(BD)
     835    else:
     836        C = 0.
     837    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
     838        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
     839        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
    819840     
    820841    BA = -BC1*SSchi+BC2*SCchi
    821842    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
    822     gam = atand2(BB,BA)
     843    gam = atan2d(BB,BA)
     844
     845    BD = (BA**2+BB**2)/rpd
     846
     847    dBAdO = 0
     848    dBAdC = -BC1*SCchi-BC2*SSchi
     849    dBAdF = BC3*SSchi
     850   
     851    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
     852    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
     853    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
     854   
     855    if BD > 0.:
     856        dGMdA = [(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD,(BA*dBBdF-BB*dBAdF)/BD]
     857    else:
     858        dGMdA = [0.0,0.0,0.0]
     859
    823860       
    824     return psi,gam
     861    return psi,gam,dPSdA,dGMdA
    825862
    826863BOH = {
    … …  
    860897}
    861898   
     899def GetKclKsl(L,N,SGLaue,psi,phi,beta):
     900    import pytexture as ptx
     901    FORPI = 12.5663706143592
     902    RSQ2PI = 0.3989422804014
     903    SQ2 = 1.414213562373
     904    Lnorm = FORPI/(2.0*L+1.)
     905    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
     906    Ksl *= RSQ2PI
     907    if SGLaue in ['m3','m3m']:
     908        Kcl = 0.0
     909        for j in range(0,L+1,4):
     910            im = j/4+1
     911            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
     912            Kcl += BOH['L='+str(l)][N-1][im-1]*pcrs*cosd(j*beta)       
     913    else:
     914        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
     915        pcrs *= RSQ2PI
     916        if N:
     917            pcrs *= SQ2
     918        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
     919            if SGLaue in ['3mR','3m1','31m']:
     920                if n%6 == 3:
     921                    Kcl = pcrs*sind(N*beta)
     922                else:
     923                    Kcl = pcrs*cosd(N*beta)
     924            else:
     925                Kcl = pcrs*cosd(N*beta)
     926        else:
     927            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
     928    return Kcl*Ksl,Lnorm
     929   
    862930def Glnh(Start,SHCoef,psi,gam,SamSym):
    863931    import pytexture as ptx
    … …  
    870938    Fln = np.zeros(len(SHCoef))
    871939    for i,term in enumerate(SHCoef):
    872          l,m,n = eval(term.strip('C'))
    873          lNorm = 4.*np.pi/(2.*l+1.)
    874          pcrs = ptx.pyplmpsi(l,m,1,psi)*RSQPI
    875          if m == 0:
    876              pcrs /= SQ2
    877          if SamSym in ['mmm',]:
    878              Ksl = pcrs*cosd(m*gam)
    879          else:
    880              Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
    881          Fln[i] = SHCoef[term]*Ksl*lNorm
     940        l,m,n = eval(term.strip('C'))
     941        lNorm = 4.*np.pi/(2.*l+1.)
     942        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
     943        pcrs *= RSQPI
     944        if m == 0:
     945            pcrs /= SQ2
     946        if SamSym in ['mmm',]:
     947            Ksl = pcrs*cosd(m*gam)
     948        else:
     949            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
     950        Fln[i] = SHCoef[term]*Ksl*lNorm
    882951    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
    883952    return ODFln
    … …  
    895964    Fln = np.zeros(len(SHCoef))
    896965    for i,term in enumerate(SHCoef):
    897          l,m,n = eval(term.strip('C'))
    898          lNorm = 4.*np.pi/(2.*l+1.)
    899          if SGData['SGLaue'] in ['m3','m3m']:
    900              Kcl = 0.0
    901              for j in range(0,l+1,4):
    902                  im = j/4+1
    903                  pcrs = ptx.pyplmpsi(l,j,1,phi)
    904                  Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
    905          else:                #all but cubic
    906              pcrs = ptx.pyplmpsi(l,n,1,phi)*RSQPI
    907              if n == 0:
    908                  pcrs /= SQ2
    909              if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
    910                 if SGData['SGLaue'] in ['3mR','3m1','31m']:
    911                     if n%6 == 3:
    912                         Kcl = pcrs*sind(n*beta)
    913                     else:
    914                         Kcl = pcrs*cosd(n*beta)
    915                 else:
    916                     Kcl = pcrs*cosd(n*beta)
    917              else:
    918                  Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
    919          Fln[i] = SHCoef[term]*Kcl*lNorm
     966        l,m,n = eval(term.strip('C'))
     967        lNorm = 4.*np.pi/(2.*l+1.)
     968        if SGData['SGLaue'] in ['m3','m3m']:
     969            Kcl = 0.0
     970            for j in range(0,l+1,4):
     971                im = j/4+1
     972                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
     973                Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
     974        else:                #all but cubic
     975            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
     976            pcrs *= RSQPI
     977            if n == 0:
     978                pcrs /= SQ2
     979            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
     980               if SGData['SGLaue'] in ['3mR','3m1','31m']:
     981                   if n%6 == 3:
     982                       Kcl = pcrs*sind(n*beta)
     983                   else:
     984                       Kcl = pcrs*cosd(n*beta)
     985               else:
     986                   Kcl = pcrs*cosd(n*beta)
     987            else:
     988                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
     989        Fln[i] = SHCoef[term]*Kcl*lNorm
    920990    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
    921991    return ODFln
    … …  
    929999        if abs(ODFln[term][1]) > 1.e-3:
    9301000            l,m,n = eval(term.strip('C'))
    931             psrs = ptx.pyplmpsi(l,m,len(psi),psi)
     1001            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
    9321002            if SamSym in ['-1','2/m']:
    9331003                if m != 0:
    … …  
    9581028                for j in range(0,l+1,4):
    9591029                    im = j/4+1
    960                     pcrs = ptx.pyplmpsi(l,j,len(beta),phi)
     1030                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
    9611031                    Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
    9621032            else:                #all but cubic
    963                 pcrs = ptx.pyplmpsi(l,n,len(beta),phi)*RSQPI
     1033                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
     1034                pcrs *= RSQPI
    9641035                if n == 0:
    9651036                    pcrs /= SQ2

  • trunk/GSASIIphsGUI.py

    r375 r379  
    19441944                    value = str(textureData['PFhkl'])
    19451945                    hkl = eval(value)
    1946                 Obj.SetValue('%d,%d,%d'%(hkl[0],hkl[1],hkl[2]))
     1946                Obj.SetValue('%d %d d'%(hkl[0],hkl[1],hkl[2]))
    19471947                textureData['PFhkl'] = hkl
    19481948            else:
    … …  
    19531953                    value = str(textureData['PFhkl'])
    19541954                    xyz = eval(value)
    1955                 Obj.SetValue('%3.1f,%3.1f,%3.1f'%(xyz[0],xyz[1],xyz[2]))
     1955                Obj.SetValue('%3.1f %3.1f %3.1f'%(xyz[0],xyz[1],xyz[2]))
    19561956                textureData['PFxyz'] = xyz
    19571957            G2plt.PlotTexture(self,data)
    … …  
    20032003            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Pole figure HKL: '),0,wx.ALIGN_CENTER_VERTICAL)
    20042004            PH = textureData['PFhkl']
    2005             pfVal = wx.TextCtrl(dataDisplay,-1,'%d,%d,%d'%(PH[0],PH[1],PH[2]),style=wx.TE_PROCESS_ENTER)
     2005            pfVal = wx.TextCtrl(dataDisplay,-1,'%d %d %d'%(PH[0],PH[1],PH[2]),style=wx.TE_PROCESS_ENTER)
    20062006        else:
    20072007            PTSizer.Add(wx.StaticText(dataDisplay,-1,' Inverse pole figure XYZ: '),0,wx.ALIGN_CENTER_VERTICAL)
    20082008            PX = textureData['PFxyz']
    2009             pfVal = wx.TextCtrl(dataDisplay,-1,'%3.1f,%3.1f,%3.1f'%(PX[0],PX[1],PX[2]),style=wx.TE_PROCESS_ENTER)
     2009            pfVal = wx.TextCtrl(dataDisplay,-1,'%3.1f %3.1f %3.1f'%(PX[0],PX[1],PX[2]),style=wx.TE_PROCESS_ENTER)
    20102010        pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue)
    20112011        pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue)
    … …  
    25342534                    mainSizer.Add(poSizer)
    25352535                    if POData[4]:
    2536                         mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
     2536                        mainSizer.Add(wx.StaticText(dataDisplay,-1,' Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
    25372537                        mainSizer.Add((0,5),0)
    25382538                        ODFSizer = wx.FlexGridSizer(2,8,2,2)

  • trunk/GSASIIstruct.py

    r378 r379  
    556556        print varstr
    557557
    558        
     558    def PrintSHPO(hapData):
     559        print '\n Spherical harmonics preferred orientation: Order:' + \
     560            str(hapData[4])+' Refine? '+str(hapData[2])
     561        ptlbls = ' names :'
     562        ptstr =  ' values:'
     563        for item in hapData[5]:
     564            ptlbls += '%12s'%(item)
     565            ptstr += '%12.4f'%(hapData[5][item])
     566        print ptlbls
     567        print ptstr
    559568   
    560569    hapDict = {}
    … …  
    603612                    hapVary.append(pfx+'MD')
    604613            else:                           #'SH' spherical harmonics
     614                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
     615                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
    605616                for item in hapData['Pref.Ori.'][5]:
    606617                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
    … …  
    639650                    Ax = hapData['Pref.Ori.'][3]
    640651                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
    641                         ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])               
     652                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
     653                else: #'SH' for spherical harmonics
     654                    PrintSHPO(hapData['Pref.Ori.'])
    642655                PrintSize(hapData['Size'])
    643656                PrintMuStrain(hapData['Mustrain'],SGData)
    … …  
    733746        print sigstr
    734747       
     748    def PrintSHPOAndSig(hapData,POsig):
     749        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
     750        ptlbls = ' names :'
     751        ptstr =  ' values:'
     752        sigstr = ' sig   :'
     753        for item in hapData[5]:
     754            ptlbls += '%12s'%(item)
     755            ptstr += '%12.4f'%(hapData[5][item])
     756            if item in POsig:
     757                sigstr += '%12.4f'%(POsig[item])
     758            else:
     759                sigstr += 12*' '
     760        print ptlbls
     761        print ptstr
     762        print sigstr
     763   
    735764    for phase in Phases:
    736765        HistoPhase = Phases[phase]['Histograms']
    … …  
    755784            else:                           #'SH' spherical harmonics
    756785                for item in hapData['Pref.Ori.'][5]:
    757                     hapData['Pref.Ori.'][5][item] = parmDict[pfx+'MD']
     786                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
    758787                    if pfx+item in sigDict:
    759788                        PhFrExtPOSig[item] = sigDict[pfx+item]
    … …  
    763792            if 'Extinction' in PhFrExtPOSig:
    764793                print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
    765             if 'MD' in PhFrExtPOSig:
    766                 print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
    767              
     794            if hapData['Pref.Ori.'][0] == 'MD':
     795                if 'MD' in PhFrExtPOSig:
     796                    print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
     797            else:
     798                PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
    768799            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
    769800                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
    … …  
    11251156    parmdict.update(zip(varylist,values))
    11261157   
    1127 def GetPrefOri(refl,G,phfx,calcControls,parmDict):
     1158def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     1159    odfCor = 1.0
     1160    H = refl[:3]
     1161    cell = G2lat.Gmat2cell(g)
     1162    Sangl = [0.,0.,0.]
     1163    if 'Bragg' in calcControls[hfx+'instType']:
     1164        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
     1165        IFCoup = True
     1166    else:
     1167        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
     1168        IFCoup = False
     1169    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1170    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
     1171    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],None,calcControls[phfx+'SHord'])
     1172    for item in SHnames:
     1173        L,N = eval(item.strip('C'))
     1174        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
     1175        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
     1176    return odfCor
     1177   
     1178def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     1179    FORPI = 12.5663706143592
     1180    odfCor = 1.0
     1181    dFdODF = {}
     1182    H = refl[:3]
     1183    cell = G2lat.Gmat2cell(g)
     1184    Sangl = [0.,0.,0.]
     1185    if 'Bragg' in calcControls[hfx+'instType']:
     1186        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
     1187        IFCoup = True
     1188    else:
     1189        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
     1190        IFCoup = False
     1191    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1192    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
     1193    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],None,calcControls[phfx+'SHord'])
     1194    for item in SHnames:
     1195        L,N = eval(item.strip('C'))
     1196        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
     1197        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
     1198        dFdODF[phfx+item] = FORPI*Kcsl
     1199    return odfCor,dFdODF
     1200   
     1201def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
    11281202    if calcControls[phfx+'poType'] == 'MD':
    11291203        MD = parmDict[phfx+'MD']
    … …  
    11361210        POcorr = sumMD/len(refl[10])
    11371211    else:   #spherical harmonics
    1138         POcorr = 1.0
     1212        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
    11391213    return POcorr
    11401214   
    1141 def GetPrefOriDerv(refl,G,phfx,calcControls,parmDict):
     1215def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
    11421216    POderv = {}
    11431217    if calcControls[phfx+'poType'] == 'MD':
    … …  
    11541228        POderv[phfx+'MD'] = sumdMD/len(refl[10])
    11551229    else:   #spherical harmonics
    1156         POcorr = 1.0
     1230        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
    11571231    return POcorr,POderv
    11581232   
    1159 def GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict):
     1233def GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
    11601234    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
    11611235    Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
    1162     Icorr *= GetPrefOri(refl,G,phfx,calcControls,parmDict)       
     1236    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)       
    11631237    return Icorr
    11641238   
    1165 def GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict):
    1166     Icorr = GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict)
     1239def GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1240    Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    11671241    pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
    1168     POcorr,dIdPO = GetPrefOriDerv(refl,G,phfx,calcControls,parmDict)
     1242    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    11691243    dIdPola /= pola
    11701244    for iPO in dIdPO:
    … …  
    13921466                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
    13931467                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
    1394                 Icorr = GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict)
     1468                Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    13951469                if 'Pawley' in Phase['General']['Type']:
    13961470                    try:
    … …  
    14881562                h,k,l = refl[:3]
    14891563                iref = pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]
    1490                 Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict)
     1564                Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    14911565                if 'Pawley' in Phase['General']['Type']:
    14921566                    try:
    … …  
    15641638def Refine(GPXfile,dlg):
    15651639    import cPickle
     1640    import pytexture as ptx
     1641    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
    15661642   
    15671643    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
    … …  
    16781754        print '\n Refinement results:'
    16791755        print 135*'-'
    1680         print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
    1681         print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
    1682         print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
     1756        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
     1757        print ' Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
     1758        print ' Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
    16831759        try:
    16841760            covMatrix = result[1]*GOF

  • trunk/fsource/pytexture.for

    r313 r379  
    44      END
    55
    6       SUBROUTINE PYPLMPSI(L,I,NPHI,PHI,PCRS)
     6      SUBROUTINE PYPLMPSI(L,I,NPHI,PHI,PCRS,DPDPS)
    77Cf2py intent(in) L
    88Cf2Py intent(in) I
    … …  
    1010Cf2py intent(in) PHI
    1111Cf2py intent(out) PCRS
    12 Cf2py depend(NPHI) PHI,PCRS
     12Cf2py intent(out) DPDPS
     13Cf2py depend(NPHI) PHI,PCRS,DPDPS
    1314
    1415      INTEGER*4     L
    … …  
    1617      REAL*4        PHI(0:NPHI-1)
    1718      REAL*4        PCRS(0:NPHI-1)
    18 
     19      REAL*4        DPDPS(0:NPHI-1)
    1920
    2021      DO K = 0,NPHI-1
    21          CALL PLMPSI(L,I,PHI(K),PCRS(K))
     22         CALL PLMPSI(L,I,PHI(K),PCRS(K),DPDPS(K))
    2223      END DO
    2324      RETURN

  • trunk/fsource/texturesubs/plmpsi.for

    r305 r379  
    1       SUBROUTINE PLMPSI(L,M,PSI,P)
     1      SUBROUTINE PLMPSI(L,M,PSI,P,DPDPS)
    22
    33!PURPOSE: Compute P(l,m,psi)
    … …  
    3131          RS = S
    3232          P = P+APR*COSD(RS*PSI)
     33          DPDPS = DPDPS-RS*APR*SIND(RS*PSI)
    3334        END DO
    3435      ELSE                                           
    … …  
    3738          RS = S
    3839          P = P+APR*SIND(RS*PSI)       
     40          DPDPS = DPDPS+RS*APR*COSD(RS*PSI)
    3941        END DO     
    4042      END IF
Note: See TracChangeset for help on using the changeset viewer.