Changeset 353


Ignore:
Timestamp:
Aug 24, 2011 4:07:29 PM (11 years ago)
Author:
vondreele
Message:

modify psvfcj.for to use sh/l only
add derivative routine to pypowder.for
GSASIIgrid.py - new LS controls for derivative type
GSASIIlattice.py - work on docs as per P Jemian
GSASIIpwd.py - major mods to use FCJpsvoight fortran code & derivatives
GSASIIpwdGUI.py - mod to use LS controls

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r346 r353  
    492492           
    493493def UpdateControls(self,data):
     494    #patch
     495    if 'deriv type' not in data:
     496        data['deriv type'] = 'analytical'
     497        data['min dM/M'] = 0.0001
     498    #end patch
    494499    '''
    495500    #Fourier controls
    … …  
    501506       
    502507    def SetStatusLine(text):
    503         Status.SetStatusText(text)
    504                                      
    505     def OnNumCycles(event):
    506         try:
    507             value = max(0,min(200,int(Ncyc.GetValue())))
    508         except ValueError:
    509             value = 3
    510         data['Ncycles'] = value
    511         Ncyc.SetValue('%d'%(value))
     508        Status.SetStatusText(text)                                     
    512509       
    513510    def OnConvergence(event):
    514511        try:
    515             value = max(0.01,min(100.,float(Cnvrg.GetValue())))
     512            value = max(1.e-9,min(1.0,float(Cnvrg.GetValue())))
    516513        except ValueError:
    517             value = 0.01
    518         data['minSumShftESD'] = value
    519         Cnvrg.SetValue('%.2f'%(value))
    520        
    521     def OnAtomShift(event):
    522         try:
    523             value = max(0.1,min(5.,float(AtShft.GetValue())))
    524         except ValueError:
    525             value = 2.0
    526         data['maxShift'] = value
    527         AtShft.SetValue('%.1f'%(value))
    528        
    529     def OnMarquardt(event):
    530         try:
    531             value = max(1.0,min(10.0,float(Marq.GetValue())))
    532         except ValueError:
    533             value = 1.0
    534         data['Marquardt'] = value
    535         Marq.SetValue('%.2f'%(value))
    536        
    537     def OnBandWidth(event):
    538         try:
    539             value = max(0,min(200,int(Band.GetValue())))
    540         except ValueError:
    541             value = 0
    542         data['bandWidth'] = value
    543         Band.SetValue('%d'%(value))
    544        
    545     def OnRestraint(event):
    546         data['restraintWeight'] = Restraint.GetValue()
    547 
     514            value = 0.0001
     515        data['min dM/M'] = value
     516        Cnvrg.SetValue('%.2g'%(value))
     517       
     518    def OnDerivType(event):
     519        data['deriv type'] = derivSel.GetValue()
     520        derivSel.SetValue(data['deriv type'])
     521       
    548522    if self.dataDisplay:
    549523        self.dataDisplay.Destroy()
    … …  
    558532    mainSizer.Add(wx.StaticText(self.dataDisplay,label=' Refinement Controls:'),0,wx.ALIGN_CENTER_VERTICAL)
    559533    LSSizer = wx.FlexGridSizer(cols=4,vgap=5,hgap=5)
    560     LSSizer.Add(wx.StaticText(self.dataDisplay,label=' Max cycles: '),0,wx.ALIGN_CENTER_VERTICAL)
    561     Ncyc = wx.TextCtrl(self.dataDisplay,-1,value='%d'%(data['Ncycles']),style=wx.TE_PROCESS_ENTER)
    562     Ncyc.Bind(wx.EVT_TEXT_ENTER,OnNumCycles)
    563     Ncyc.Bind(wx.EVT_KILL_FOCUS,OnNumCycles)
    564     LSSizer.Add(Ncyc,0,wx.ALIGN_CENTER_VERTICAL)
    565     LSSizer.Add(wx.StaticText(self.dataDisplay,label=' Min sum(shift/esd)^2: '),0,wx.ALIGN_CENTER_VERTICAL)
    566     Cnvrg = wx.TextCtrl(self.dataDisplay,-1,value='%.2f'%(data['minSumShftESD']),style=wx.TE_PROCESS_ENTER)
     534    LSSizer.Add(wx.StaticText(self.dataDisplay,label='Refinement derivatives: '),0,wx.ALIGN_CENTER_VERTICAL)
     535    Choice=['analytic','numeric']
     536    derivSel = wx.ComboBox(parent=self.dataDisplay,value=data['deriv type'],choices=Choice,
     537        style=wx.CB_READONLY|wx.CB_DROPDOWN)
     538    derivSel.SetValue(data['deriv type'])
     539    derivSel.Bind(wx.EVT_COMBOBOX, OnDerivType)   
     540    LSSizer.Add(derivSel,0,wx.ALIGN_CENTER_VERTICAL)
     541    LSSizer.Add(wx.StaticText(self.dataDisplay,label=' Min delta-M/M: '),0,wx.ALIGN_CENTER_VERTICAL)
     542    Cnvrg = wx.TextCtrl(self.dataDisplay,-1,value='%.2g'%(data['min dM/M']),style=wx.TE_PROCESS_ENTER)
    567543    Cnvrg.Bind(wx.EVT_TEXT_ENTER,OnConvergence)
    568544    Cnvrg.Bind(wx.EVT_KILL_FOCUS,OnConvergence)
    569545    LSSizer.Add(Cnvrg,0,wx.ALIGN_CENTER_VERTICAL)
    570     LSSizer.Add(wx.StaticText(self.dataDisplay,label=' Max atom shift: '),0,wx.ALIGN_CENTER_VERTICAL)
    571     AtShft = wx.TextCtrl(self.dataDisplay,-1,value='%.1f'%(data['maxShift']),style=wx.TE_PROCESS_ENTER)
    572     AtShft.Bind(wx.EVT_TEXT_ENTER,OnAtomShift)
    573     AtShft.Bind(wx.EVT_KILL_FOCUS,OnAtomShift)
    574     LSSizer.Add(AtShft,0,wx.ALIGN_CENTER_VERTICAL)
    575     LSSizer.Add(wx.StaticText(self.dataDisplay,label=' Marquardt factor: '),0,wx.ALIGN_CENTER_VERTICAL)
    576     Marq = wx.TextCtrl(self.dataDisplay,-1,value='%.2f'%(data['Marquardt']),style=wx.TE_PROCESS_ENTER)
    577     Marq.Bind(wx.EVT_TEXT_ENTER,OnMarquardt)
    578     Marq.Bind(wx.EVT_KILL_FOCUS,OnMarquardt)
    579     LSSizer.Add(Marq,0,wx.ALIGN_CENTER_VERTICAL)
    580     LSSizer.Add(wx.StaticText(self.dataDisplay,label=' Matrix band width: '),0,wx.ALIGN_CENTER_VERTICAL)
    581     Band = wx.TextCtrl(self.dataDisplay,-1,value='%d'%(data['bandWidth']),style=wx.TE_PROCESS_ENTER)
    582     Band.Bind(wx.EVT_TEXT_ENTER,OnBandWidth)
    583     Band.Bind(wx.EVT_KILL_FOCUS,OnBandWidth)
    584     LSSizer.Add(Band,0,wx.ALIGN_CENTER_VERTICAL)
    585     Restraint = wx.CheckBox(self.dataDisplay,-1,label='Modify restraint weights?')
    586     Restraint.Bind(wx.EVT_CHECKBOX, OnRestraint)
    587     Restraint.SetValue(data['restraintWeight'])
    588     LSSizer.Add(Restraint,0,wx.ALIGN_CENTER_VERTICAL)
    589546   
    590547   
    … …  
    753710                data = {
    754711                    #least squares controls
    755                     'Ncycles':3,'maxShift':2.0,'bandWidth':0,'Marquardt':1.0,'restraintWeight':False,
    756                     'minSumShftESD':0.01,
     712                    'deriv type':'analytic','min dM/M':0.0001,
    757713                    #Fourier controls
    758714                    'mapType':'Fobs','d-max':100.,'d-min':0.2,'histograms':[],

  • trunk/GSASIIlattice.py

    r342 r353  
    2323
    2424def sec2HMS(sec):
     25    """Convert time in sec to H:M:S string
     26   
     27    :param sec: time in seconds
     28    return: H:M:S string (to nearest 100th second)
     29   
     30    """
    2531    H = int(sec/3600)
    2632    M = int(sec/60-H*60)
    … …  
    2935   
    3036def rotdMat(angle,axis=0):
    31     '''Prepare rotation matrix for angle in degrees about axis(=0,1,2)
    32     Returns numpy 3,3 array
    33     '''
     37    """Prepare rotation matrix for angle in degrees about axis(=0,1,2)
     38
     39    :param angle: angle in degrees
     40    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
     41    :return: rotation matrix - 3x3 numpy array
     42
     43    """
    3444    if axis == 2:
    3545        return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
    … …  
    4050       
    4151def rotdMat4(angle,axis=0):
    42     '''Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
    43     Returns numpy 4,4 array
    44     '''
     52    """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
     53
     54    :param angle: angle in degrees
     55    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
     56    :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling)
     57
     58    """
    4559    Mat = rotdMat(angle,axis)
    4660    return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
    4761   
    4862def fillgmat(cell):
    49     '''Compute lattice metric tensor from unit cell constants
    50     cell is tuple with a,b,c,alpha, beta, gamma (degrees)
    51     returns 3x3 numpy array
    52     '''
     63    """Compute lattice metric tensor from unit cell constants
     64
     65    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
     66    :return: 3x3 numpy array
     67
     68    """
    5369    a,b,c,alp,bet,gam = cell
    5470    g = np.array([
    … …  
    5975           
    6076def cell2Gmat(cell):
    61     '''Compute real and reciprocal lattice metric tensor from unit cell constants
    62     cell is tuple with a,b,c,alpha, beta, gamma (degrees)
    63     returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
    64     '''
     77    """Compute real and reciprocal lattice metric tensor from unit cell constants
     78
     79    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
     80    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
     81
     82    """
    6583    g = fillgmat(cell)
    6684    G = nl.inv(g)       
    … …  
    6886
    6987def A2Gmat(A):
    70     '''Fill reciprocal metric tensor (G) from A
    71     returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
    72     '''
     88    """Fill real & reciprocal metric tensor (G) from A
     89
     90    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
     91    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
     92
     93    """
    7394    G = np.zeros(shape=(3,3))
    7495    G = [
    … …  
    80101
    81102def Gmat2A(G):
    82     'Extract A from reciprocal metric tensor (G)'
     103    """Extract A from reciprocal metric tensor (G)
     104
     105    :param G: reciprocal maetric tensor (3x3 numpy array
     106    :return: A = [G11,G22,G33,2*G12,2*G13,2*G23]
     107
     108    """
    83109    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
    84110   
    85111def cell2A(cell):
     112    """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters
     113
     114    :param cell: [a,b,c,alpha,beta,gamma] (degrees)
     115    :return: G reciprocal metric tensor as 3x3 numpy array
     116
     117    """
    86118    G,g = cell2Gmat(cell)
    87119    return Gmat2A(G)
    88120
    89121def A2cell(A):
    90     '''Compute unit cell constants from A tensor
    91     returns tuple with a,b,c,alpha, beta, gamma (degrees)
    92     '''
     122    """Compute unit cell constants from A
     123
     124    :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor
     125    :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters
     126
     127    """
    93128    G,g = A2Gmat(A)
    94129    return Gmat2cell(g)
    95130
    96131def Gmat2cell(g):
    97     '''Compute lattice parameters from real metric tensor (g)
    98     returns tuple with a,b,c,alpha, beta, gamma (degrees)
    99     Alternatively,compute reciprocal lattice parameters from inverse metric tensor (G)
    100     returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
    101     '''
     132    """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G)
     133    The math works the same either way.
     134
     135    :param g (or G): real (or reciprocal) metric tensor 3x3 array
     136    :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees)
     137
     138    """
    102139    oldset = np.seterr('raise')
    103140    a = np.sqrt(max(0,g[0][0]))
    … …  
    111148
    112149def invcell2Gmat(invcell):
    113     '''Compute real and reciprocal lattice metric tensor from reciprocal
     150    """Compute real and reciprocal lattice metric tensor from reciprocal
    114151       unit cell constants
    115     invcell is tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
    116     returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
    117     '''
     152       
     153    :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees)
     154    :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
     155
     156    """
    118157    G = fillgmat(invcell)
    119158    g = nl.inv(G)
    … …  
    121160       
    122161def calc_rVsq(A):
    123     'Compute the square of the reciprocal lattice volume (V* **2) from A'
     162    """Compute the square of the reciprocal lattice volume (V* **2) from A'
     163
     164    """
    124165    G,g = A2Gmat(A)
    125166    rVsq = nl.det(G)
    … …  
    129170   
    130171def calc_rV(A):
    131     'Compute the reciprocal lattice volume (V*) from A'
     172    """Compute the reciprocal lattice volume (V*) from A
     173    """
    132174    return np.sqrt(calc_rVsq(A))
    133175   
    134176def calc_V(A):
    135     'Compute the real lattice volume (V) from A'
     177    """Compute the real lattice volume (V) from A
     178    """
    136179    return 1./calc_rV(A)
    137180
    138181def A2invcell(A):
    139     '''Compute reciprocal unit cell constants from A
     182    """Compute reciprocal unit cell constants from A
    140183    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
    141     '''
     184    """
    142185    G,g = A2Gmat(A)
    143186    return Gmat2cell(G)
    144187
    145188def cell2AB(cell):
    146     '''Computes orthogonalization matrix from unit cell constants
     189    """Computes orthogonalization matrix from unit cell constants
    147190    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
    148191    returns tuple of two 3x3 numpy arrays (A,B)
    149192       A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
    150193       B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B*x) = x
    151     '''
     194    """
    152195    G,g = cell2Gmat(cell)
    153196    cellstar = Gmat2cell(G)
    … …  
    164207   
    165208def U6toUij(U6):
    166     '''Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
     209    """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
    167210    returns
    168     '''
     211    """
    169212    U = np.zeros(shape=(3,3))
    170213    U = [
    … …  
    175218       
    176219def Uij2betaij(Uij,G):
    177     '''
     220    """
    178221    Convert Uij to beta-ij tensors
    179222    input:
    … …  
    182225    returns:
    183226    beta-ij - numpy array [beta-ij]
    184     '''
     227    """
    185228    pass
    186229   
    187230def CosSinAngle(U,V,G):
    188     ''' calculate sin & cos of angle betwee U & V in generalized coordinates
     231    """ calculate sin & cos of angle betwee U & V in generalized coordinates
    189232    defined by metric tensor G
    190233    input:
    … …  
    193236    return:
    194237        cos(phi) & sin(phi)
    195     '''
     238    """
    196239    u = U/nl.norm(U)
    197240    v = V/nl.norm(V)
    … …  
    201244   
    202245def criticalEllipse(prob):
    203     '''
     246    """
    204247    Calculate critical values for probability ellipsoids from probability
    205     '''
     248    """
    206249    if not ( 0.01 <= prob < 1.0):
    207250        return 1.54
    … …  
    212255   
    213256def CellBlock(nCells):
    214     '''
     257    """
    215258    Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1
    216259    currently only works for nCells = 0 or 1 (not >1)
    217     '''
     260    """
    218261    if nCells:
    219262        N = 2*nCells+1
    … …  
    241284#   
    242285def _combinators(_handle, items, n):
    243     ''' factored-out common structure of all following combinators '''
     286    """ factored-out common structure of all following combinators """
    244287    if n==0:
    245288        yield [ ]
    … …  
    250293            yield this_one + cc
    251294def combinations(items, n):
    252     ''' take n distinct items, order matters '''
     295    """ take n distinct items, order matters """
    253296    def skipIthItem(items, i):
    254297        return items[:i] + items[i+1:]
    255298    return _combinators(skipIthItem, items, n)
    256299def uniqueCombinations(items, n):
    257     ''' take n distinct items, order is irrelevant '''
     300    """ take n distinct items, order is irrelevant """
    258301    def afterIthItem(items, i):
    259302        return items[i+1:]
    260303    return _combinators(afterIthItem, items, n)
    261304def selections(items, n):
    262     ''' take n (not necessarily distinct) items, order matters '''
     305    """ take n (not necessarily distinct) items, order matters """
    263306    def keepAllItems(items, i):
    264307        return items
    265308    return _combinators(keepAllItems, items, n)
    266309def permutations(items):
    267     ''' take all items, order matters '''
     310    """ take all items, order matters """
    268311    return combinations(items, len(items))
    269312
    … …  
    361404                                   
    362405def GetBraviasNum(center,system):
    363     '''Determine the Bravais lattice number, as used in GenHBravais
    364          center = one of: P, C, I, F, R (see SGLatt from GSASIIspc.SpcGroup)
    365          lattice = is cubic, hexagonal, tetragonal, orthorhombic, trigonal (R)
    366              monoclinic, triclinic (see SGSys from GSASIIspc.SpcGroup)
    367        Returns a number between 0 and 13
    368           or throws an exception if the setting is non-standard
    369        '''
     406    """Determine the Bravais lattice number, as used in GenHBravais
     407   
     408    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
     409    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
     410             'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
     411    :return: a number between 0 and 13
     412          or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
     413    """
    370414    if center.upper() == 'F' and system.lower() == 'cubic':
    371415        return 0
    … …  
    399443
    400444def GenHBravais(dmin,Bravais,A):
    401     '''Generate the positionally unique powder diffraction reflections
    402     input:
    403        dmin is minimum d-space
    404        Bravais is 0-13 to indicate lattice type (see GetBraviasNum)
    405        A is reciprocal cell tensor (see Gmat2A or cell2A)
    406     returns:
    407        a list of tuples containing: h,k,l,d-space,-1   
    408     '''
    409 # Bravais in range(14) to indicate Bravais lattice:
    410 #   0 F cubic
    411 #   1 I cubic
    412 #   2 P cubic
    413 #   3 R hexagonal (trigonal not rhombohedral)
    414 #   4 P hexagonal
    415 #   5 I tetragonal
    416 #   6 P tetragonal
    417 #   7 F orthorhombic
    418 #   8 I orthorhombic
    419 #   9 C orthorhombic
    420 #  10 P orthorhombic
    421 #  11 C monoclinic
    422 #  12 P monoclinic
    423 #  13 P triclinic
    424 # A - as defined in calc_rDsq
    425 # returns HKL = [h,k,l,d,0] sorted so d largest first
     445    """Generate the positionally unique powder diffraction reflections
     446     
     447    :param dmin: minimum d-spacing in A
     448    :param Bravais: lattice type (see GetBraviasNum)
     449        Bravais is one of::
     450             0 F cubic
     451             1 I cubic
     452             2 P cubic
     453             3 R hexagonal (trigonal not rhombohedral)
     454             4 P hexagonal
     455             5 I tetragonal
     456             6 P tetragonal
     457             7 F orthorhombic
     458             8 I orthorhombic
     459             9 C orthorhombic
     460            10 P orthorhombic
     461            11 C monoclinic
     462            12 P monoclinic
     463            13 P triclinic
     464    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
     465    :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
     466           
     467    """
    426468    import math
    427469    if Bravais in [9,11]:
    … …  
    512554   
    513555def GenHLaue(dmin,SGData,A):
    514     '''Generate the crystallographically unique powder diffraction reflections
     556    """Generate the crystallographically unique powder diffraction reflections
    515557    for a lattice and Bravais type
    516     Input:
    517         dmin - minimum d-spacing
    518         SGData - space group dictionary with at least:
    519             SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
    520                      '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
    521             SGLatt - lattice centering = 'P','A','B','C','I','F'
    522             SGUniq - code for unique monoclinic axis = 'a','b','c'
    523         A - 6 terms as defined in calc_rDsq
    524     Return;
    525         HKL = list of [h,k,l,d] sorted with largest d first and is unique
     558   
     559    :param dmin: minimum d-spacing
     560    :param SGData: space group dictionary with at least::
     561   
     562        'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
     563                 '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
     564        'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
     565        'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m')
     566            otherwise ' '
     567       
     568    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
     569    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
    526570            part of reciprocal space ignoring anomalous dispersion
    527     '''
     571           
     572    """
    528573    import math
    529574    SGLaue = SGData['SGLaue']

  • trunk/GSASIIpwd.py

    r350 r353  
    492492fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    493493               
    494 def getBackground(pfx,parmDict,bakType,xdata):
    495     yb = np.zeros_like(xdata)
    496     if bakType == 'chebyschev':
    497         iBak = 0
    498         while True:
    499             key = pfx+'Back:'+str(iBak)
    500             try:
    501                 yb += parmDict[key]*(xdata-xdata[0])**iBak
    502                 iBak += 1
    503             except KeyError:
    504                 break
    505     return yb
    506 
    507494def getWidths(pos,sig,gam,shl):
    508495    widths = [np.sqrt(sig)/100.,gam/200.]
    … …  
    535522    return intens*Df(xdata)*DX/dx
    536523
     524def getBackground(pfx,parmDict,bakType,xdata):
     525    yb = np.zeros_like(xdata)
     526    if bakType == 'chebyschev':
     527        iBak = 0
     528        while True:
     529            key = pfx+'Back:'+str(iBak)
     530            try:
     531                yb += parmDict[key]*(xdata-xdata[0])**iBak
     532                iBak += 1
     533            except KeyError:
     534                break
     535    return yb
     536   
     537def getBackgroundDerv(pfx,parmDict,bakType,xdata):
     538    dydb = []
     539    if bakType == 'chebyschev':
     540        iBak = 0
     541        while True:
     542            if pfx+'Back:'+str(iBak) in parmDict:
     543                dydb.append((xdata-xdata[0])**iBak)
     544                iBak += 1
     545            else:
     546                break
     547    return dydb
     548
    537549#use old fortran routine
    538 def getFCJVoigt3(pos,intens,sig,gam,shl,xdata):
     550def getFCJVoigt3(pos,sig,gam,shl,xdata):
    539551   
    540552    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
    541553    Df /= np.sum(Df)
    542     return intens*Df
     554    return Df
     555
     556def getdFCJVoigt3(pos,sig,gam,shl,xdata):
     557   
     558    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) #might have to make these numpy arrays?
     559    sumDf = np.sum(Df)
     560    return Df/sumDf,dFdp,dFds,dFdg,dFdsh
     561   
    543562
    544563def getPeakProfile(parmDict,xdata,varyList,bakType):
    … …  
    552571    X = parmDict['X']
    553572    Y = parmDict['Y']
    554     shl = max(parmDict['SH/L'],0.0005)
     573    shl = max(parmDict['SH/L'],0.002)
    555574    Ka2 = False
    556575    if 'Lam1' in parmDict.keys():
    … …  
    589608            elif not iBeg-iFin:     #peak above high limit
    590609                return yb+yc
    591             yc[iBeg:iFin] += getFCJVoigt3(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
     610            yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
    592611            if Ka2:
    593612                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
    … …  
    596615                iFin = min(lenX,iFin+kdelt)
    597616                if iBeg-iFin:
    598                     yc[iBeg:iFin] += getFCJVoigt3(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
     617                    yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
    599618            iPeak += 1
    600619        except KeyError:        #no more peaks to process
    601620            return yb+yc
    602 
     621           
     622def getPeakProfileDerv(parmDict,xdata,varyList,bakType):
     623# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
     624    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
     625    if 'Back:0' in varyList:            #background derivs are in front if present
     626        dMdb = getBackgroundDerv('',parmDict,bakType,xdata)
     627        dMdv[0:len(dMdb)] = dMdb
     628       
     629    dx = xdata[1]-xdata[0]
     630    U = parmDict['U']
     631    V = parmDict['V']
     632    W = parmDict['W']
     633    X = parmDict['X']
     634    Y = parmDict['Y']
     635    shl = max(parmDict['SH/L'],0.002)
     636    Ka2 = False
     637    if 'Lam1' in parmDict.keys():
     638        Ka2 = True
     639        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
     640        kRatio = parmDict['I(L2)/I(L1)']
     641    iPeak = 0
     642    while True:
     643        try:
     644            pos = parmDict['pos'+str(iPeak)]
     645            intens = parmDict['int'+str(iPeak)]
     646            sigName = 'sig'+str(iPeak)
     647            tanth = tand(pos/2.0)
     648            costh = cosd(pos/2.0)
     649            if sigName in varyList:
     650                sig = parmDict[sigName]
     651            else:
     652                sig = U*tanth**2+V*tanth+W
     653                dsdU = tanth**2
     654                dsdV = tanth
     655                dsdW = 1.0
     656            sig = max(sig,0.001)          #avoid neg sigma
     657            gamName = 'gam'+str(iPeak)
     658            if gamName in varyList:
     659                gam = parmDict[gamName]
     660            else:
     661                gam = X/costh+Y*tanth
     662                dgdX = 1.0/costh
     663                dgdY = tanth
     664            gam = max(gam,0.001)             #avoid neg gamma
     665            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
     666            iBeg = np.searchsorted(xdata,pos-fmin)
     667            lenX = len(xdata)
     668            if not iBeg:
     669                iFin = np.searchsorted(xdata,pos+fmin)
     670            elif iBeg == lenX:
     671                iFin = iBeg
     672            else:
     673                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
     674            if not iBeg+iFin:       #peak below low limit
     675                iPeak += 1
     676                continue
     677            elif not iBeg-iFin:     #peak above high limit
     678                break
     679            dMdpk = np.zeros(shape=(6,len(xdata)))
     680            dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
     681            for i in range(5):
     682                dMdpk[i][iBeg:iFin] = 100.*dx*intens*dMdipk[i]
     683            dMdpk[0][iBeg:iFin] = dMdipk[0]
     684            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
     685            if Ka2:
     686                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     687                kdelt = int((pos2-pos)/dx)               
     688                iBeg = min(lenX,iBeg+kdelt)
     689                iFin = min(lenX,iFin+kdelt)
     690                if iBeg-iFin:
     691                    dMdipk = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
     692                    for i in range(5):
     693                        dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk[i]
     694                    dMdpk[0][iBeg:iFin] += kRatio*dMdipk[0]
     695                    dMdpk[5][iBeg:iFin] += dMdipk[0]
     696                    dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
     697            for parmName in ['pos','int','sig','gam']:
     698                try:
     699                    idx = varyList.index(parmName+str(iPeak))
     700                    dMdv[idx] = dervDict[parmName]
     701                except ValueError:
     702                    pass
     703            if 'U' in varyList:
     704                dMdv[varyList.index('U')] += dsdU*dervDict['sig']
     705            if 'V' in varyList:
     706                dMdv[varyList.index('V')] += dsdV*dervDict['sig']
     707            if 'W' in varyList:
     708                dMdv[varyList.index('W')] += dsdW*dervDict['sig']
     709            if 'X' in varyList:
     710                dMdv[varyList.index('X')] += dgdX*dervDict['gam']
     711            if 'Y' in varyList:
     712                dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
     713            if 'SH/L' in varyList:
     714                dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
     715            if 'I(L2)/I(L1)' in varyList:
     716                dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
     717            iPeak += 1
     718        except KeyError:        #no more peaks to process
     719            break
     720    return dMdv
     721       
    603722def Dict2Values(parmdict, varylist):
    604723    '''Use before call to leastsq to setup list of values for the parameters
    … …  
    611730    parmdict.update(zip(varylist,values))
    612731   
    613 def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
     732def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False,controls=None):
    614733   
    615734    def SetBackgroundParms(Background):
    … …  
    651770        insVary = []
    652771        for i,flag in enumerate(insFlags):
    653             if flag:
     772            if flag and insNames[i] in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
    654773                insVary.append(insNames[i])
    655774        instDict = dict(zip(insNames,insVals))
    656775        instDict['X'] = max(instDict['X'],0.01)
    657776        instDict['Y'] = max(instDict['Y'],0.01)
    658         instDict['SH/L'] = max(instDict['SH/L'],0.0001)
     777        instDict['SH/L'] = max(instDict['SH/L'],0.002)
    659778        return dataType,instDict,insVary
    660779       
    … …  
    743862                    ptstr += 12*' '
    744863            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
     864               
     865    def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):
     866        parmdict.update(zip(varylist,values))
     867        return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType)
    745868           
    746869    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
    747870        parmdict.update(zip(varylist,values))
    748         M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
     871        M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)-ydata)
    749872        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
    750873        if dlg:
    … …  
    754877        return M
    755878       
    756     Ftol = 0.0001
     879    if controls:
     880        Ftol = controls['min dM/M']
     881        derivType = controls['deriv type']
     882    else:
     883        Ftol = 0.0001
     884        derivType = 'analytic'
    757885    if oneCycle:
    758886        Ftol = 1.0
    … …  
    778906            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    779907            try:
    780                 result = so.leastsq(errPeakProfile,values,full_output=True,epsfcn=1.e-8,ftol=Ftol,
    781                     args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     908                if derivType == 'analytic':
     909                    result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
     910                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     911                    ncyc = int(result[2]['nfev']/2)
     912                else:
     913                    result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,
     914                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     915                    ncyc = int(result[2]['nfev']/len(varyList))
    782916            finally:
    783917                dlg.Destroy()
    784918            runtime = time.time()-begin   
    785919            chisq = np.sum(result[2]['fvec']**2)
    786             ncyc = int(result[2]['nfev']/len(varyList))
    787920            Values2Dict(parmDict, varyList, result[0])
    788921            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
    … …  
    9291062        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
    9301063        }
     1064    global parmDict2
     1065    parmDict2 = {
     1066        'pos0':5.7,'int0':10.0,'sig0':0.5,'gam0':1.0,
     1067        'U':1.,'V':-1.,'W':0.1,'X':0.0,'Y':2.,'SH/L':0.004,
     1068        'Back0':5.,'Back1':-0.02,'Back2':.004,
     1069        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
     1070        }
    9311071    global varyList
    9321072    varyList = []
    … …  
    9491089    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
    9501090   
     1091def test2(name,delt):
     1092    if NeedTestData: TestData()
     1093    varyList = [name,]
     1094    xdata = np.linspace(5.6,5.8,800)
     1095    hplot = plotter.add('derivatives test for '+name).gca()
     1096    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
     1097    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
     1098    parmDict2[name] += delt
     1099    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
     1100    hplot.plot(xdata,(y1-y0)/delt,'r+')
     1101   
    9511102   
    9521103
    … …  
    9551106    global plotter
    9561107    plotter = plot.PlotNotebook()
    957     test0()
     1108#    test0()
     1109    test2('I(L2)/I(L1)',.0001)
    9581110    print "OK"
    9591111    plotter.StartEventLoop()

  • trunk/GSASIIpwdGUI.py

    r345 r353  
    9696    def OnPeakFit(FitPgm,oneCycle=False):
    9797        SaveState()
    98         print 'Peak Fitting:'
     98        controls = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,self.root, 'Controls'))
     99        if not controls:
     100            controls = {'deriv type':'analytic','min dM/M':0.0001,}     #fill in defaults if needed
     101        print 'Peak Fitting with '+controls['deriv type']+' derivatives:'
    99102        PatternId = self.PatternId
    100103        PickId = self.PickId
    … …  
    109112        wx.BeginBusyCursor()
    110113        try:
    111             G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data,oneCycle)
     114            G2pwd.DoPeakFit(FitPgm,peaks,background,limits,inst,data,oneCycle,controls)
    112115        finally:
    113116            wx.EndBusyCursor()   

  • trunk/fsource/powsubs/psvfcj.for

    r210 r353  
    1       SUBROUTINE PSVFCJ(DTT,TTHETA,SL,HL,SIG,GAM,PRFUNC,DPRDT,SLPART,
    2      1  HLPART,SIGPART,GAMPART)
     1      SUBROUTINE PSVFCJ(DTT,TTHETA,SIG,GAM,SHL,PRFUNC,
     2     1  DPRDT,SIGPART,GAMPART,SHLPART)
    33
    44!PURPOSE: Compute function & derivatives for Pseudovoigt profile
    … …  
    88!   [L.W. Finger, D.E. Cox & A.P. Jephcoat (1994) J. Appl. Cryst.,27,892-900.]
    99! coded 11/95 by B. H. Toby (NIST). revised version
    10 ! parameterized as asym1=S/L asym2=H/L
     10! parameterized as asym=S/L+H/L
    1111
    1212
    … …  
    1717      REAL*4        DTT                 !delta 2-theta in centidegrees                 
    1818      REAL*4        TTHETA              !2-theta in centidegrees             
    19       REAL*4        SL,HL               !S/L & H/L               
    2019      REAL*4        SIG,GAM             
     20      REAL*4        SHL                 !S/L + H/L               
    2121      REAL*4        PRFUNC             
    2222      REAL*4        DPRDT               
    23       REAL*4        SLPART,HLPART       
    2423      REAL*4        SIGPART,GAMPART     
     24      REAL*4        SHLPART     
    2525
    2626!INCLUDE STATEMENTS:
    … …  
    3737      REAL*4        DFDA               
    3838      REAL*4        DGDA               
    39       REAL*4        DGDB               
    4039      REAL*4        DYDA               
    41       REAL*4        DYDB               
    4240      REAL*4        SIN2THETA2          ! sin(2theta)**2
    4341      REAL*4        COS2THETA           ! cos(2theta)
    … …  
    4947      REAL*4        COSDELTA2           ! cos(Delta)**2
    5048      REAL*4        A                   ! asym1 [coff(7)]
    51       REAL*4        B                   ! asym2 [coff(8)]
    5249      REAL*4        APB                 ! (A+B)
    53       REAL*4        AMB                 ! (A-B)
    5450      REAL*4        APB2                ! (A+B)**2
    5551      REAL*4        TTHETAD             ! Two Theta in degrees
    … …  
    6460      REAL*4        SUMWDGDA            !      sum of w dGdA
    6561      REAL*4        SUMWRDGDA           !      sum of w R dGdA
    66       REAL*4        SUMWDGDB            !      sum of w dGdB
    67       REAL*4        SUMWRDGDB           !      sum of w R dGdB
    6862      REAL*4        SUMWGDRD2T          !      sum of w G dRd(2theta)
    6963      REAL*4        SUMWGDRDSIG         !      sum of w G dRdp(n)
    7064      REAL*4        SUMWGDRDGAM         !      sum of w G dRdp(n)
    7165      REAL*4        SUMWGDRDA           
    72       REAL*4        SUMWGDRDB           
    7366      REAL*4        EMIN                ! 2phi minimum
    7467      REAL*4        EINFL               ! 2phi of inflection point
    … …  
    126119C Asymmetry terms
    127120c
    128       A = SL            ! A = S/L in FCJ
    129       B = HL            ! B = H/L in FCJ
    130       ApB = A+B
    131       AmB = A-B
     121      A = SHL/2.            ! A = S/L in FCJ
     122      ApB = SHL
    132123      ApB2 = ApB*ApB
    133124c
    134125C handle the case where there is asymmetry
    135126c
    136       IF (A .ne. 0.0 .or. B .ne. 0.0) then
    137         Einfl = Acosd(SQRT(1.0 + AmB**2)*cos2THETA) ! 2phi(infl) FCJ eq 5 (degrees)
     127      IF (A .ne. 0.0) then
     128        Einfl = Acosd(cos2THETA) ! 2phi(infl) FCJ eq 5 (degrees)
    138129        tmp2 = 1.0 + ApB2
    139130        tmp = SQRT(tmp2)*cos2THETA
    140131c
    141 C Treat case where A or B is zero - Set Einfl = 2theta
    142 c
    143         if (A.eq.0.0 .or. B .eq. 0.0)Einfl = Acosd(cos2THETA)
     132C Treat case where A is zero - Set Einfl = 2theta
     133c
     134        if (A.eq.0.0) Einfl = Acosd(cos2THETA)
    144135        if (abs(tmp) .le. 1.0) then
    145136          Emin = Acosd(tmp)      ! 2phi(min) FCJ eq 4 (degrees)
    … …  
    154145        endif
    155146        if (tmp1 .gt. 0 .and. abs(tmp) .le. 1.0) then
    156           dEmindA = -ApB*cos2THETA/SQRT(tmp1) ! N. B. symm w/r A,B
     147          dEmindA = -ApB*cos2THETA/SQRT(tmp1)
    157148        ELSE
    158149          dEmindA = 0.0
    … …  
    203194        sumWdGdA = 0.
    204195        sumWRdGdA = 0.
    205         sumWdGdB = 0.
    206         sumWRdGdB = 0.
    207196        sumWGdRd2t = 0.
    208197        sumWGdRdsig = 0.
    … …  
    241230c       
    242231          if(abs(delta-emin) .gt. abs(einfl-emin))then
    243             if ( A.ge.B) then
    244232c
    245233C N.B. this is the only place where d()/dA <> d()/dB
    246234c
    247               G = 2.0*B*F*RcosDELTA
    248               dGdA = 2.0*B*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
    249               dGdB = dGdA + 2.0*F*RcosDELTA
    250             else
    251235              G = 2.0*A*F*RcosDELTA
    252               dGdB = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
    253               dGdA = dGdB + 2.0*F*RcosDELTA
    254             endif
     236              dGdA = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
    255237          else                                            ! delta .le. einfl .or. min(A,B) .eq. 0
    256238            G = (-1.0 + ApB*F) * RcosDELTA
    257239            dGdA = RcosDELTA*(F - tanDELTA*dDELTAdA
    258240     1             + ApB*F*tanDELTA*dDELTAdA + ApB*dFdA)
    259             dGdB = dGdA
    260241          endif
    261242
    … …  
    265246          sumWdGdA = sumWdGdA + wp(k+it) * dGdA
    266247          sumWRdGdA = sumWRdGdA + wp(k+it) * R * dGdA
    267           sumWdGdB = sumWdGdB + wp(k+it) * dGdB
    268           sumWRdGdB = sumWRdGdB + wp(k+it) * R * dGdB
    269248          sumWGdRd2t = sumWGdRd2t + WG * dRdT ! N.B. 1/centidegrees
    270249          sumWGdRdsig = sumWGdRdsig + WG * dRdS
    … …  
    278257        dydA = (-(sumWRG*sumWdGdA) +
    279258     1    sumWG*(sumWRdGdA- 100.0 * todeg * sumWGdRdA)) * RsumWG2
    280         dydB = (-(sumWRG*sumWdGdB) +
    281      1    sumWG*(sumWRdGdB- 100.0 * todeg * sumWGdRdA)) * RsumWG2
    282259        sigpart = sumWGdRdsig * RsumWG
    283260        gampart = sumWGdRdgam * RsumWG
    … …  
    289266        CALL PSVOIGT(DTT,SIG,GAM,R,dRdT,dRdS,dRdG)
    290267        PRFUNC = R
    291         dydA = 0.002 * sign(1.0,TTHETA - DTT)
    292         dydB = dydA
     268        dydA = 0.004 * sign(1.0,TTHETA - DTT)
    293269        sigpart = dRdS
    294270        gampart = dRdG
    295271        DPRDT = -dRdT
    296272      END IF
    297       SLPART = DYDA
    298       HLPART = DYDB
     273      SHLPART = DYDA
    299274     
    300275      RETURN

  • trunk/fsource/pypowder.for

    r352 r353  
    33C TTHETA in degrees
    44C SPH is S/L + H/L
     5C RETURNS FUNCTION ONLY
    56Cf2py intent(in) NPTS
    67Cf2py intent(in) DTT
    … …  
    1415
    1516      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
    16       SL = SPH/2.0
    1717      FW = (2.355*SQRT(SIG)+GAM)/100.0
    1818      FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
    1919      FMAX = 15.0*FW
    2020      DO I=0,NPTS-1
    21         CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SL,SL,SIG,GAM,
    22      1    PRFUNC(I),DPRDT,SLPART,HLPART,SIGPART,GAMPART)
     21        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
     22     1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART)
    2323      END DO
    2424      RETURN
    2525      END
     26
     27      SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC,
     28     1  DPRDT,SIGPART,GAMPART,SLPART)
     29C DTT in degrees
     30C TTHETA in degrees
     31C SPH is S/L + H/L
     32C RETURNS FUNCTION & DERIVATIVES
     33Cf2py intent(in) NPTS
     34Cf2py intent(in) DTT
     35cf2py depend(NPTS) DTT
     36Cf2py intent(in) TTHETA
     37Cf2py intent(in) SIG
     38Cf2py intent(in) GAM
     39Cf2py intent(in) SPH
     40Cf2py intent(out) PRFUNC
     41Cf2py depend(NPTS) PRFUNC
     42Cf2py intent(out) DPRDT
     43Cf2py depend(NPTS) DPRDT
     44Cf2py intent(out) SIGPART
     45Cf2py depend(NPTS) SIGPART
     46Cf2py intent(out) GAMPART
     47Cf2py depend(NPTS) GAMPART
     48Cf2py intent(out) SLPART
     49Cf2py depend(NPTS) SLPART
     50
     51      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
     52     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
     53      FW = (2.355*SQRT(SIG)+GAM)/100.0
     54      FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
     55      FMAX = 15.0*FW
     56      DO I=0,NPTS-1
     57        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
     58     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I))
     59        DPRDT(I) = DPRDT(I)*100.
     60      END DO
     61      RETURN
     62      END
     63
Note: See TracChangeset for help on using the changeset viewer.