Changeset 379
- Timestamp:
- Sep 22, 2011 10:11:18 AM (14 years ago)
- Location:
- trunk
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified trunk/GSASIIlattice.py ¶
r378 r379 788 788 789 789 def 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. 790 802 if IFCoup: 791 803 GSomeg = sind(Gangls[2]+Tth) … … 817 829 BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg 818 830 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)] 819 840 820 841 BA = -BC1*SSchi+BC2*SCchi 821 842 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 823 860 824 return psi,gam 861 return psi,gam,dPSdA,dGMdA 825 862 826 863 BOH = { … … 860 897 } 861 898 899 def 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 862 930 def Glnh(Start,SHCoef,psi,gam,SamSym): 863 931 import pytexture as ptx … … 870 938 Fln = np.zeros(len(SHCoef)) 871 939 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 882 951 ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln)))) 883 952 return ODFln … … 895 964 Fln = np.zeros(len(SHCoef)) 896 965 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 920 990 ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln)))) 921 991 return ODFln … … 929 999 if abs(ODFln[term][1]) > 1.e-3: 930 1000 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) 932 1002 if SamSym in ['-1','2/m']: 933 1003 if m != 0: … … 958 1028 for j in range(0,l+1,4): 959 1029 im = j/4+1 960 pcrs = ptx.pyplmpsi(l,j,len(beta),phi)1030 pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi) 961 1031 Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta) 962 1032 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 964 1035 if n == 0: 965 1036 pcrs /= SQ2 -
TabularUnified trunk/GSASIIphsGUI.py ¶
r375 r379 1944 1944 value = str(textureData['PFhkl']) 1945 1945 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])) 1947 1947 textureData['PFhkl'] = hkl 1948 1948 else: … … 1953 1953 value = str(textureData['PFhkl']) 1954 1954 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])) 1956 1956 textureData['PFxyz'] = xyz 1957 1957 G2plt.PlotTexture(self,data) … … 2003 2003 PTSizer.Add(wx.StaticText(dataDisplay,-1,' Pole figure HKL: '),0,wx.ALIGN_CENTER_VERTICAL) 2004 2004 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) 2006 2006 else: 2007 2007 PTSizer.Add(wx.StaticText(dataDisplay,-1,' Inverse pole figure XYZ: '),0,wx.ALIGN_CENTER_VERTICAL) 2008 2008 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) 2010 2010 pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue) 2011 2011 pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue) … … 2534 2534 mainSizer.Add(poSizer) 2535 2535 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) 2537 2537 mainSizer.Add((0,5),0) 2538 2538 ODFSizer = wx.FlexGridSizer(2,8,2,2) -
TabularUnified trunk/GSASIIstruct.py ¶
r378 r379 556 556 print varstr 557 557 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 559 568 560 569 hapDict = {} … … 603 612 hapVary.append(pfx+'MD') 604 613 else: #'SH' spherical harmonics 614 controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4] 615 controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5]) 605 616 for item in hapData['Pref.Ori.'][5]: 606 617 hapDict[pfx+item] = hapData['Pref.Ori.'][5][item] … … 639 650 Ax = hapData['Pref.Ori.'][3] 640 651 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.']) 642 655 PrintSize(hapData['Size']) 643 656 PrintMuStrain(hapData['Mustrain'],SGData) … … 733 746 print sigstr 734 747 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 735 764 for phase in Phases: 736 765 HistoPhase = Phases[phase]['Histograms'] … … 755 784 else: #'SH' spherical harmonics 756 785 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] 758 787 if pfx+item in sigDict: 759 788 PhFrExtPOSig[item] = sigDict[pfx+item] … … 763 792 if 'Extinction' in PhFrExtPOSig: 764 793 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) 768 799 SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]], 769 800 'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]], … … 1125 1156 parmdict.update(zip(varylist,values)) 1126 1157 1127 def GetPrefOri(refl,G,phfx,calcControls,parmDict): 1158 def 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 1178 def 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 1201 def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): 1128 1202 if calcControls[phfx+'poType'] == 'MD': 1129 1203 MD = parmDict[phfx+'MD'] … … 1136 1210 POcorr = sumMD/len(refl[10]) 1137 1211 else: #spherical harmonics 1138 POcorr = 1.01212 POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict) 1139 1213 return POcorr 1140 1214 1141 def GetPrefOriDerv(refl,G, phfx,calcControls,parmDict):1215 def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): 1142 1216 POderv = {} 1143 1217 if calcControls[phfx+'poType'] == 'MD': … … 1154 1228 POderv[phfx+'MD'] = sumdMD/len(refl[10]) 1155 1229 else: #spherical harmonics 1156 POcorr = 1.01230 POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict) 1157 1231 return POcorr,POderv 1158 1232 1159 def GetIntensityCorr(refl,G, phfx,hfx,calcControls,parmDict):1233 def GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): 1160 1234 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity 1161 1235 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) 1163 1237 return Icorr 1164 1238 1165 def GetIntensityDerv(refl,G, phfx,hfx,calcControls,parmDict):1166 Icorr = GetIntensityCorr(refl,G, phfx,hfx,calcControls,parmDict)1239 def GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): 1240 Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) 1167 1241 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) 1169 1243 dIdPola /= pola 1170 1244 for iPO in dIdPO: … … 1392 1466 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift 1393 1467 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) 1395 1469 if 'Pawley' in Phase['General']['Type']: 1396 1470 try: … … 1488 1562 h,k,l = refl[:3] 1489 1563 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) 1491 1565 if 'Pawley' in Phase['General']['Type']: 1492 1566 try: … … 1564 1638 def Refine(GPXfile,dlg): 1565 1639 import cPickle 1640 import pytexture as ptx 1641 ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics 1566 1642 1567 1643 def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): … … 1678 1754 print '\n Refinement results:' 1679 1755 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) 1683 1759 try: 1684 1760 covMatrix = result[1]*GOF -
TabularUnified trunk/fsource/pytexture.for ¶
r313 r379 4 4 END 5 5 6 SUBROUTINE PYPLMPSI(L,I,NPHI,PHI,PCRS )6 SUBROUTINE PYPLMPSI(L,I,NPHI,PHI,PCRS,DPDPS) 7 7 Cf2py intent(in) L 8 8 Cf2Py intent(in) I … … 10 10 Cf2py intent(in) PHI 11 11 Cf2py intent(out) PCRS 12 Cf2py depend(NPHI) PHI,PCRS 12 Cf2py intent(out) DPDPS 13 Cf2py depend(NPHI) PHI,PCRS,DPDPS 13 14 14 15 INTEGER*4 L … … 16 17 REAL*4 PHI(0:NPHI-1) 17 18 REAL*4 PCRS(0:NPHI-1) 18 19 REAL*4 DPDPS(0:NPHI-1) 19 20 20 21 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)) 22 23 END DO 23 24 RETURN -
TabularUnified trunk/fsource/texturesubs/plmpsi.for ¶
r305 r379 1 SUBROUTINE PLMPSI(L,M,PSI,P )1 SUBROUTINE PLMPSI(L,M,PSI,P,DPDPS) 2 2 3 3 !PURPOSE: Compute P(l,m,psi) … … 31 31 RS = S 32 32 P = P+APR*COSD(RS*PSI) 33 DPDPS = DPDPS-RS*APR*SIND(RS*PSI) 33 34 END DO 34 35 ELSE … … 37 38 RS = S 38 39 P = P+APR*SIND(RS*PSI) 40 DPDPS = DPDPS+RS*APR*COSD(RS*PSI) 39 41 END DO 40 42 END IF
Note: See TracChangeset
for help on using the changeset viewer.