Changeset ac995be in sasmodels


Ignore:
Timestamp:
Jun 12, 2018 10:29:02 AM (4 years ago)
Author:
Adam Washington <adam.washington@…>
Branches:
master
Children:
b032200
Parents:
d7af1c6
Message:

Get log sesans roughly working

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/sesans.py

    rd7af1c6 rac995be  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi  # type: ignore 
    16 from scipy.special import j0, j1 
     16from scipy.special import j1 
    1717 
    1818 
     
    3939 
    4040    # transform arrays 
    41     _H = None  # type: np.ndarray 
    42     _H0 = None # type: np.ndarray 
     41    _H = None   # type: np.ndarray 
     42    _H0 = None  # type: np.ndarray 
    4343 
    4444    def __init__(self, z, SElength, lam, zaccept, Rmax): 
    4545        # type: (np.ndarray, float, float) -> None 
    46         #import logging; logging.info("creating SESANS transform") 
    4746        self.q = z 
    4847        self._set_hankel(SElength, lam, zaccept, Rmax) 
     
    5251        G0 = np.dot(self._H0, Iq) 
    5352        G = np.dot(self._H.T, Iq) 
     53        G0 = G[0]  # FIXME: This is a kludge 
    5454        P = G - G0 
    5555        return P 
     
    5757    def _set_hankel(self, SElength, lam, zaccept, Rmax): 
    5858        # type: (np.ndarray, float, float) -> None 
    59         # Force float32 arrays, otherwise run into memory problems on some machines 
     59        # Force float32 arrays, otherwise run into memory problems on 
     60        # some machines 
    6061        SElength = np.asarray(SElength, dtype='float32') 
    6162 
     
    6465        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 
    6566        # q = np.arange(q_min, q_max, q_min, dtype='float32') 
    66         q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2), 
    67                              dtype=np.float32)) 
    68         q = np.hstack([0], q) 
     67        # q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2), 
     68        #                      dtype=np.float32)) 
     69        q = np.exp(np.linspace(np.log(q_min), np.log(q_max), 10*SElength.size, 
     70                               dtype=np.float32)) 
     71        q = np.hstack([[0], q]) 
    6972 
    70         H0 = np.pi * (q[1:]**2 - q[-1]**2) 
     73        H0 = np.pi * (q[1:]**2 - q[:-1]**2) * (q[1:] - q[:-1]) 
    7174 
    7275        # repq = np.tile(q, (SElength.size, 1)).T 
    7376        H = np.outer(q, SElength) 
    7477        j1(H, out=H) 
    75         H *= q 
     78        H *= q.reshape((-1, 1)) 
    7679        H = H[1:] - H[:-1] 
    7780        H *= 2 * np.pi / SElength 
    7881 
    7982        lam = np.asarray(lam, dtype=np.float32) 
    80         reptheta = np.outer(q, lam) 
     83        reptheta = np.outer(q[1:], lam) 
    8184        reptheta /= np.float32(2*np.pi) 
    8285        np.arcsin(reptheta, out=reptheta) 
    8386        # reptheta = np.arcsin(repq*replam/2*np.pi) 
    8487        mask = reptheta > zaccept 
    85         H[mask] = 0 
     88        # H[mask] = 0 
    8689 
    8790        # H = np.zeros((q.size, SElength.size), dtype=np.float32) 
    8891        # H0 = q * 0 
    89         assert(H.shape == (q.size, SElength.size)) 
     92        assert(H.shape == (q.size-1, SElength.size)) 
    9093 
    91         self.q_calc = q 
     94        self.q_calc = q[1:] 
    9295        self._H, self._H0 = H, H0 
Note: See TracChangeset for help on using the changeset viewer.