new libraries (devel)

#!/usr/bin/env python

import numpy as np

from pysen.config import coverage

from pysen.messung import make_messung

if __name__ == "__main__":

    qtau_inp ='messung_18286-SampleA_9A.inp'

    params = {

        'temperature_controller'      : 'julabo_biooven',

        'temperature_controller_ref'  : 'julabo_biooven',

    }

    qscan_def = (

        #mode        Q     L      count  detlim  atten  revolver

        {

            'mode': 'standard',

            'qmin': 0.10,

            'lmax': 8e-10,

            'counts': 450,

            'det_angle_limit': 5.0,

            'attenuator': 'open',

            'revolver'  : 3,

            'taulist'   : np.linspace(1,11,6),

        },

        {

            'mode': 'standard',

            'qmin': 0.14,

            'lmax': 11e-10,

            'counts': 450,

            'det_angle_limit': 5.0,

            'attenuator': 'open',

            'revolver'  : 1,

            'taulist'   : np.linspace(1,11,6),

        },

    )

    fd = open(qtau_inp)

    glbs = globals()

    exec(compile(fd.read(), qtau_inp, 'exec'), glbs)

    print(len(data))

    for d in data:

        for k in d:

            if k=='limits': continue

            if k=='plot_data': continue

            print(k, d[k])

    #make_messung("messung_test", *qscan_def, **params)

#!/usr/bin/env python

from __future__ import print_function

import sys

import os.path

import warnings

import numpy as np

import matplotlib.pyplot as plt

from   matplotlib import cm

from   scipy.optimize import curve_fit, OptimizeWarning

import libstapler as ls

def main(*args, **kwargs):

    model = kwargs.pop('model', 'KWW-Norm') # fit model (see fit_sqt)

    mpars = kwargs.pop('model_pars', {} )   # model parameters

    norm  = kwargs.pop('norm' , False  )    # whether to 'scale' results or not

    rings = kwargs.pop('rings', (0,))       # which q-rings to use

    taus  = kwargs.pop('taus', None)        # tau range

    plot_summary = kwargs.pop('plot_summary', False)

    data_dir     = kwargs.pop('data_dir', '.')

    results = []           # list of fit results

    ncolor  = 3

    print("# ", model)

    #for sample in args:

    for sample in args:

        ncolor = max(len(sample[1:]), 3)

        label  = sample[0]

        print("# ", label)

        if model.startswith('Diffusion'):

            print("# q     dq      D          dD         chi2  extra_pars")

        elif model.startswith('Zilman-Granek'):

            print("# q     dq      Gamma      dGamma     chi2  extra_pars")

        else:

            print("# q     dq      tau0       dtau0      chi2  extra_pars")

        plt.figure(figsize=(8,8))

        fitres = []

        for norm_fac, files, slabel in sample[1:]:

            files = files.split()

            # nice filename

            outfile    = label.replace('/','_').replace(' ','_').replace('=','').lower()

            # 1st figure S(q,t) versus q

            for data_index in rings:

                q, dq, tsam, data = stitch(*files, data_index=data_index, data_dir=data_dir,

                                           norm_fac=norm_fac*int(norm))

                tau  = data[0]

                if taus:

                    sel  = np.logical_and(taus[0]<=tau, tau<=taus[1])

                    data = data[:,sel]

                tau  = data[0]

                sqt  = abs(data[1])

                serr = abs(data[2])

                if not len(tau):

                    continue

                #alabel = r"Q=(%.3f$\pm$%.3f) $\AA^{-1}$" % (q, dq)

                alabel = r"Q=%.3f $\AA^{-1}$" % q

                if slabel:

                    alabel = r"%s %s" % (alabel, slabel)

                else:

                    alabel = r"%s T=%.0fK" % (alabel, tsam)

                plt.errorbar(tau, sqt, yerr=serr, fmt='s', markersize=5, label=alabel) # data

                try:

                    par0, epar0, fit_func, chi2, popt = fit_sqt(q, data, model=model, **mpars)

                    ftau = np.logspace(np.log10(min(tau))-0.2, np.log10(max(tau))+0.2)

                    fsqt = fit_func(ftau)

                    plt.plot(ftau, fsqt, '--', color=acolor) # fit

                    fitres.append((q, dq, par0, epar0, tsam))

                    print("%.4f, %.4f, %.3e, %.3e, %.3g,\t" % (q, dq, par0, epar0, chi2), end="")

                    if len(popt)>1:

                        print(", ".join(["%.3g" % _x for _x in popt[1:]]), end=",")

                    print(tsam)

                except (TypeError, IndexError) as e:

                    print("data fitting error: %s " % files)

                except (RuntimeError, OptimizeWarning) as e:

                    print("data fitting error: %s (%s)" % (files, e))

                    plt.plot(tau, sqt, '--', color=acolor) # data

                if slabel:

                    flabel = slabel.replace('/','_').replace(' ','_').replace('=','').lower()

                    write_csv("%s-q-%.3f-%s.csv" % (outfile,q,flabel), data, label=label, q=q, dq=dq, temp=tsam)

                else:

                    write_csv("%s-q-%.3f.csv" % (outfile,q), data, label=label, q=q, dq=dq, temp=tsam)

        results.append((label, fitres))

        #plt.xscale('log')

        #plt.yscale('log')

        #plt.xlim(xmin= 0.01, xmax=100.0)

        plt.ylim(ymin= 0.5, ymax=1.1)

        plt.legend(loc='best')

        plt.xlabel(r'$\tau$ [ns]')

        plt.ylabel(r'$S(Q,\tau)/S(Q,0)$')

        plt.title("%s (%s)" % (label, model))

        plt.grid(which='both')

        plt.savefig("%s-sqt.pdf" % str(outfile))

    if plot_summary and model != 'None':

        plt.figure(figsize=(8,8))

        color  = get_color(len(results))

        for label, fitres in results:

            if len(fitres)<1: continue

            fitres = np.asarray(fitres).T

            q     = fitres[0]

            dq    = fitres[1]

            par0  = fitres[2]

            dpar0 = fitres[3]

            tsam  = fitres[4]

            acolor = next(color)

            #plt.errorbar(q, par0, yerr=dpar0, xerr=dq,

            #             fmt='s', label=label, color=acolor)

            #plt.xlabel(r'$Q \; (\AA^{-1})$')

            plt.errorbar(tsam, par0, yerr=dpar0, xerr=1.0,

                         fmt='s', label=label, color=acolor)

            plt.xlabel(r'T [K]')

            ## par0 = A*q^-2

            #Q0 = 2.5

            #pa0 = np.average(np.log(par0)+2*np.log(q))

            #A = np.exp(pa0)

            #B = -2

            #print(1/A,B)

            #xq = np.linspace(0.8*min(q), 1.2*Q0)

            #plt.plot(xq, A*xq**B, '--', color=acolor,

            #         label=r'Q$^{%.3g}$/%.3g' % (B,1/A))

        plt.title("Model: %s" % model)

        #plt.yscale('log')

        #plt.xlim(xmin=0.0)

        #plt.ylim(ymin=0.0)

        if model.lower().startswith('diffusion'):

            plt.ylabel(r'D$\ \; (\AA^{2} ns^{-1}) $')

        elif model.lower().startswith('zilman'):

            plt.ylabel(r'$\Gamma \; (ns^{-1}) $')

        else:

            plt.ylabel(r'$\tau_0 \ \; (ns) $')

        plt.legend(loc='best')

        plt.grid()

        plt.savefig(model.lower()+'-summary.pdf')

    plt.show()

if __name__ == "__main__":

    sample1_fs = [ r'Sample1 FirstScan',

                # q=0.10, T=310->295

                (1.00, "b__7840t__4", ""),

                (1.00, "b__7841t__4 b__7844t__4", ""),

                (1.00, "b__7842t__4", ""),

                (1.00, "b__7843t__4", ""),

                # q=0.40 T=295K

                (1.00, "b__7846t__4", ""),

                # q=0.56 T=295K

                (1.00, "b__7845t__4", ""),

            ]

    sample1_t390 = [ r'Sample1 T390K',

                #(1.00, "b__7851t__4", ""), #Q=0.14

                #(1.00, "b__7850t__4", ""), #Q=0.20

                (1.00, "b__7849t__4", ""), #Q=0.28

                (1.00, "b__7848t__4", ""), #Q=0.40

            ]

    sample1_t327 = [ r'Sample1 T327K',

                (1.00, "b__7867t__4", ""),             #Q=0.20

                (1.00, "b__7866t__4 b__7870t__4", ""), #Q=0.40

                (1.00, "b__7869t__4", ""),             #Q=0.56

               ]

    sample1_q04 = [ r'Sample1 Q=0.47',

                (1.00, "b__7848t__4", ""),             #Q=0.4, T=390K

                #(1.00, "b__7866t__4 b__7870t__4", ""), #Q=0.4, T=327K

                (1.025, "b__7939t__4 ", "T=305K"),     #Q=0.4, T=305K

                (1.025, "b__7922t__4 ", "T=310K"),     #Q=0.4, T=310K

               ]

    sample1_q04_res = [ r'Sample1 Q=0.5 (New Resolution)',

                #(1.000, "b__7939t__4        ", "T=305K/TiZr"),    #Q=0.4, T=305K

                #(1.000, "b__7922t__4        ", "T=310K/TiZr"),    #Q=0.4, T=310K

                #(1.000, "b__7938t__4        ", "T=315K/TiZr"),    #Q=0.4, T=315K

                (0.98, "b__7939t__4.newres ", "T=305K/10K"),     #Q=0.4, T=305K

                (0.98, "b__7922t__4.newres ", "T=310K/10K"),     #Q=0.4, T=310K

                (0.98, "b__7938t__4.newres ", "T=315K/10K"),     #Q=0.4, T=315K

                ]

    betaKWW = np.power(0.49,1/1.23)

    main(

        #sample1_fs,

        #sample1_t390,

        #sample1_t327,

        #sample1_q04,

        sample1_q04_res,

        #

        #model='KWW',

        #model_pars=dict(bounds=([0.0,0.00],[1e4,5.0])),

        model='KWW-Fix',

        model_pars={'beta': betaKWW},

        #model_pars={'beta': 2.7},

        #norm=True,

        plot_summary=True,

        rings=(0,),

        data_dir='nobgr',

        )

#!/usr/bin/env python

"""

Stapler library module

"""

from __future__ import print_function

import os.path

import ast

from collections import OrderedDict

import numpy as np

from scipy.optimize import curve_fit

#import warnings

#from scipy.optimize import OptimizeWarning

#warnings.filterwarnings('ignore', '', RuntimeWarning)

#warnings.filterwarnings('error' , '', OptimizeWarning)

A2ns = 1e-7 # A^2/ns = 1e-7 cm^2/s

#function_library = {

#    'KWW-Norm'       : (3, lambda t, t0, beta, A:    A*np.exp(-(t/t0)**beta)), # b-fitting

#    'KWW'            : (2, lambda t, t0, beta:         np.exp(-(t/t0)**beta)), # b-fitting

#    'KWW-Fix'        : (1, lambda t, t0:               np.exp(-(t/t0)**beta)), # np b-fitting

#    'KWW-Fix+Const'  : (2, lambda t, t0, A:      (1-A)*np.exp(-(t/t0)**beta)+A), #

#    'Zilman-Granek'  : (1, lambda t, G:                np.exp(-(G*t)**beta) ), # no b-fitting

#    'Exp'            : (1, lambda t, t0:               np.exp(-(t/t0))      ), #

#    'Exp-Norm'       : (2, lambda t, t0, A:          A*np.exp(-(t/t0))      ), #

#    'Exp+Const'      : (2, lambda t, t0, A:      (1-A)*np.exp(-(t/t0))+A    ), #

#    'Diffusion'      : (1, lambda t, G:                np.exp(-G*q*q*t)     ), # simple diffusion

#    'Diffusion-Norm' : (1, lambda t, G,  A:          A*np.exp(-G*q*q*t)     ), # simple diffusion

#    'Diffusion+Const': (2, lambda x, G,  A: (1-A)*np.exp(-G*q*q*x)+A        ), # simple diffusion

#    'Power'          : (2, lambda x, b,  A:     A*x**b               ), #

#    'Linear'         : (2, lambda x, x0, A:     A*x-x0               ), #

#    'Constant'       : (1, lambda x, A:         A*np.ones_like(x)    ), #S

#    'Cumulant'       : (2, lambda x, c1, c2:    np.exp(-c1*x+c2*x**2/2.0)   ),

#    'Diffusion+2ND'  : (3, lambda t, G,  c2:    np.exp(-G*q*q*t+c2*t*t) ), # diffusion

#}

def kww(tau, tau0, beta):

    "Kolrausch-Williams-Watson"

    return np.exp(-(tau/tau0)**beta)

def diffusion(tau, DQ2):

    "Simple Diffusion"

    return np.exp(-DQ2*tau)

def zilman_granek(tau, gamma):

    "Zilman-Granek"

    beta = 2/3

    return np.exp(-(gamma*tau)**beta)

def fit_sqt(data, model=zilman_granek, **kwargs):

    """

    fit S(Q,t)/S(Q,0) according to a model

    returns a tuple:

        (parameter, parameter_error, fit_function, chisq, pars)

    """

    #A       = kwargs.pop('A',  1.0)

    #G       = kwargs.pop('G',  1.0)

    #t0      = kwargs.pop('t0', 1.0)

    #beta    = kwargs.pop('beta', 2.0/3.0)

    bounds  = kwargs.pop('bounds', (-np.inf, np.inf))

    max_tau = kwargs.pop('max_tau', np.inf)

    tau  = data[0]

    sqt  = data[1]

    serr = data[2]

    popt, pcov = curve_fit(model, tau[tau<max_tau], sqt[tau<max_tau], sigma=serr[tau<max_tau],

                           bounds=bounds, maxfev=10000, ftol=1e-5)

    par, epar = popt[0], np.sqrt(pcov[0,0])

    fit_func = lambda x: model(x, *popt)

    chisq = np.sum(((model(tau[tau<max_tau]) - sqt[tau<max_tau])/serr[tau<max_tau])**2)

    chisq = chisq/len(sqt[tau<max_tau] - len(par))

    return (par,epar, fit_func, chisq, popt)

def delta_q(q, meta,  nt=42):

    """

    Calculate delta(q)

    """

    try:

        l    = meta['lambda']

        lmin = meta['lam_min']

        lmax = meta['lam_max']

        t1   = meta['tbin_lower']

        t2   = meta['tbin_upper']

        sint2 = 0.25*q*l/np.pi        # sin(theta/2.0)

        dlam  = np.abs(lmax-lmin)     # lambda band (max)

        l1  = lmax - (nt-t1)*dlam/nt  # selected lambda band (lower)

        l2  = lmax - (nt-t2)*dlam/nt  # selected lambda band (upper)

        dq1 = np.abs(4.0*np.pi/l1*sint2 - q)

        dq2 = np.abs(4.0*np.pi/l2*sint2 - q)

        return np.sqrt(dq1**2 + dq2**2)/2.0

    except KeyError:

        return np.zeros_like(q)

def stitch(*files, **keyw):

    """

    Stitch data from several files

    """

    data_dir    = keyw.pop('data_dir', '.')

    data_index  = keyw.pop('data_index', 0)

    norm_fac    = keyw.pop('norm_fac',   0)

    data = np.zeros((3,0))  #nr columns in a b-file

    qa   = []

    sq   = []

    ts   = []

    for filename in files:

        filename = os.path.join(data_dir, filename)

        if not os.path.exists(filename):

            continue

        m, d = read_wbfile(filename)[data_index]

        try:

            d = d[:3,:]

            q = m['q']

            t = m.get('temp_act') or m.get('temp')

            ts.append(t)

            qa.append(q)

            sq.append(delta_q(q, m))

            data = np.concatenate((data,d),axis=1)

        except IndexError:

            print("data reading error", filename)

    qa = np.asarray(qa)

    sq = np.asarray(sq)

    ts = np.asarray(ts)

    if norm_fac:

        data[1] = data[1]/norm_fac

        data[2] = data[2]/norm_fac

    data = data[:, data[0].argsort()]

    return np.average(qa), np.sqrt(np.sum(sq**2)), np.average(ts), data

def _evaluate(value):

    "evaluate value"

    try:

        value = ast.literal_eval(value) # convert to numbers

    except (ValueError, SyntaxError):

        pass

    return value

def read_wbfile(filename):

    """

    Read b/w files produced by echodet

    Returns list of (metadata, data) pairs

    """

    result = []             # result

    data   = []             # data

    meta   = OrderedDict()  # meta data

    isdata = False

    iline  = 0

    with open(filename, 'rt') as fd:

        for line in fd.readlines():

            line = line.strip()

            if not line:

                continue

            iline += 1

            if iline <= 2:  # first two lines are info

                meta['info%d' % iline] = line

                continue

            if line  == 'values':  # values indicate start of data block

                isdata = True

                continue

            if line == '#eod' or line == '#nxt':

                result.append((meta,np.asarray(data).T))

                data = []

                meta = OrderedDict()

                isdata = False

                iline  = 0

                if line == '#eod':

                    break

                continue

            xl = line.split()

            if isdata:

                try:

                    data.append([float(x) for x in xl])

                except ValueError:

                    pass

            else:

                key   = xl[0]

                value = " ".join(xl[1:])

                meta[key] = _evaluate(value)

    return result

def write_csv(filename, data, **kwargs):

    "write processed csv file"

    header = ["%s" % kwargs.pop('label'),]

    for key in kwargs:

        header.append("%s, %s" % (key, kwargs.get(key)))

    header.append(" ")

    header.append("%13.13s%13.13s%13.13s" % ("tau","S(Q,t)","dS(Q,t)"))

    np.savetxt(filename, data.T, fmt='%13.6f',header="\n".join(header), delimiter=',')

"""

==================================

NSE scattering evaluator for a Rouse polymer chain

==================================

-----------------------------------------------------------------------

M. Monkenbusch

Juelich Center of Neutron Science

Forschungszentrum Juelich

m.monkenbusch@fz-juelich.de

-----------------------------------------------------------------------

This is a little toolbox to estimate the effect of contrast and

H/D content on the dynamics data of a polymer melt in the Rouse

regime as obtained by NSE spectroscopy.

It is not intended to scrutinize Rouse behaviour and deviations from that,

for that the time functions should be computed with the methods

given e.g. in the book of Doi-Edwards.

But it will give a first impression on scattering contributions etc.

In particular:

consideration concerning the roles of coherent an incoherent scattering

in polymer samples

please observe:

-----------------------------------------------------------

units for

Q                                        --> 1/Angstroem

lambda                               --> Angstroem

scattering length                 --> cm

scattering length density    --> cm**-2

scattering cross section      --> cm**2

neutron flux                       --> neutrons/cm**2/s

molecular weights              --> g/mol

density                              --> g/cm**3

sample dimensions:

thichkness,d                      --> cm

area                                   --> cm**2

Rouse rate parameter

Wl4                                    --> Angstroem**4/nanoseconds

t                                        --> nanoseconds

------------------------------------------------------------

"""

import numpy as np

from numpy import pi, exp, sqrt

from scipy.special import erfc

from scipy.integrate import quad

_eps  = 1e-33

_fac  = 2/sqrt(pi)

def debye(x):

    "Debye function f(x)"

    return 2*(exp(-x)+1+x)/x**2

# Rouse kernel $h(u) = \frac{2}{\pi} \int_0^\infty \cos(ux)/x^2 (1-e^{-x^2}) dx$

#              $h(u) = \frac{2}{\sqrt{pi}} \exp{-(u/2)^2} - u erfc(u/2) $

def h_rouse(u):

    "Rouse kernel h(u)"

    u2 = u/2

    return _fac*exp(-(u2)**2) - u*erfc(u2)

# exp(-(u+sqrt(omegat)*h_rouse(u/sqrt(omegat))))

def skernel_rouse(u, xot):

    "Rouse S(Q,t) inegrand"

    #xot  = omegat #sqrt(omegat)

    if xot<_eps:

        return exp(-u)  # h(oo) = 0

    return exp(-(u+xot*h_rouse(u/xot)))

#  $ S(Q,t)/S(Q,0) = \int_0^\infty dx exp(-u + \sqrt(omegat)

def sqt_rouse(omegat):

    "Rouse dynamic structure factor"

    return np.asarray([ quad(skernel_rouse, 0, np.inf, epsrel=1e-6,args=(_ot,)) for _ot in omegat ])