[group] [eFit] no-fit yes-simulate color=blue name=muScatt-Eq1 no-log
HF
[name]
muscatt-example-5t1
[number parameters]
17,
[description]
MuScatt@QtiSAS
MuScatt: Multiple Scattering Program for Small Angle Scattering
[Henrich Frielinghaus]
https://jugit.fz-juelich.de/sans/muscatt
Integration of MuScatt c-code to Fitting interface of QtiSAS
(Fit-Compile & Fit-Table)
[Vitaliy Pipic]
http://iffwww.iff.kfa-juelich.de/~pipich/dokuwiki/lib/exe/fetch.php/qtisas/download/muscatt_qtisas.pdf
v. 2019-10-29
[x]
Q
[y]
I
[parameter names]
lowQ,bgksub,lambda,d,transm,sigmia,hgqmod,hgqcut,resmod,dlfwhm,detlen,sampap,collen,entrap,reswdt,thkmod,activeOutput,
[initial values]
0,0,5,0.1,0.554,3.6442,0,0.0689,-1,0.1,20,0.01,20,0.03,2,0,1
[adjustibility]
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0
[parameter description]
switch for lowQ extrapolation: 0 automatic 1 guinier 2 power law 3 constant,incoherent background subtraction: -1 none 0 automatic (real)value>0 gives exact level to be subtracted,lambda (Angstroem) ,sample thickness in cm ,experimental transmission as measured,sigma due to incoherent scattering and absorption calculated from table (will be divided by 4pi for full 3d angle),high q extrapolation (0 for off 1 for on),q range from which pure power law is achieved,resolution mode #n (-1 is off;; 0 for no combination);; combines with DEMUX/MUX calculation of type #n (1-6),fraction of lambda resolution at fwhm;; i.e. delta lambda / lambda,detector distance [m] for lowest Q measurements,sample aperture [m],collimation distance [m] for lowest Q measurements;; negative values if not specified (symmetric setting),entrance aperture [m];; negative values if not specified (symmetric setting) ,resolution cutoff (in terms of multiples of sigma),thickness correction mode #n (<=0 is off;; combines with DEMUX/MUX calculation of type #n (1-6) and (7-12) for (1-6) with resolution correction,0..8: 0 input data;; 1 Eq.1;; 2 Eq. 3;; 3 Eq. 2...,
[h-headers]
#include
#include
#include
#include "IncludedFunctions/muscatt-alloc.h"
#include "IncludedFunctions/muscatt-fftsg.h"
#include "IncludedFunctions/muscatt-functions.h"
[included functions]
static bool calculated=false;
static double In[300][16];
static double Qc[300];
static double III[300];
static double dIII[300];
static int orderQ[300];
//+++
//+++ saveToTable
//+++
void saveToTable(int nd, void * ParaM )
{
//+++ col. names
( (struct functionT *) ParaM)->tableColNames= (char**)malloc(19 * sizeof(char *));
((struct functionT *) ParaM)->tableColNames[0]="Q";
((struct functionT *) ParaM)->tableColNames[1]="I";
((struct functionT *) ParaM)->tableColNames[2]="dI";
((struct functionT *) ParaM)->tableColNames[3]="Eq1";
((struct functionT *) ParaM)->tableColNames[4]="dEq1";
((struct functionT *) ParaM)->tableColNames[5]="Eq3";
((struct functionT *) ParaM)->tableColNames[6]="dEq3";
((struct functionT *) ParaM)->tableColNames[7]="Eq2";
((struct functionT *) ParaM)->tableColNames[8]="dEq2";
((struct functionT *) ParaM)->tableColNames[9]="Eq1reversal";
((struct functionT *) ParaM)->tableColNames[10]="dEq1reversal";
((struct functionT *) ParaM)->tableColNames[11]="Eq3reversal";
((struct functionT *) ParaM)->tableColNames[12]="dEq3reversal";
((struct functionT *) ParaM)->tableColNames[13]="Eq2reversal";
((struct functionT *) ParaM)->tableColNames[14]="dEq2reversal";
((struct functionT *) ParaM)->tableColNames[15]="CorrIntWithReso";
((struct functionT *) ParaM)->tableColNames[16]="dCorrIntWithReso";
((struct functionT *) ParaM)->tableColNames[17]="CorrIntWithThickness";
((struct functionT *) ParaM)->tableColNames[18]="dCorrIntWithThickness";
// column designations (X, Y, etc.)
((struct functionT *) ParaM)->tableColDestinations= new int[19];
((struct functionT *) ParaM)->tableColDestinations[0]=1;
((struct functionT *) ParaM)->tableColDestinations[1]=2;
((struct functionT *) ParaM)->tableColDestinations[2]=5;
((struct functionT *) ParaM)->tableColDestinations[3]=2;
((struct functionT *) ParaM)->tableColDestinations[4]=5;
((struct functionT *) ParaM)->tableColDestinations[5]=2;
((struct functionT *) ParaM)->tableColDestinations[6]=5;
((struct functionT *) ParaM)->tableColDestinations[7]=2;
((struct functionT *) ParaM)->tableColDestinations[8]=5;
((struct functionT *) ParaM)->tableColDestinations[9]=2;
((struct functionT *) ParaM)->tableColDestinations[10]=5;
((struct functionT *) ParaM)->tableColDestinations[11]=2;
((struct functionT *) ParaM)->tableColDestinations[12]=5;
((struct functionT *) ParaM)->tableColDestinations[13]=2;
((struct functionT *) ParaM)->tableColDestinations[14]=5;
((struct functionT *) ParaM)->tableColDestinations[15]=2;
((struct functionT *) ParaM)->tableColDestinations[16]=5;
((struct functionT *) ParaM)->tableColDestinations[17]=2;
((struct functionT *) ParaM)->tableColDestinations[18]=5;
((struct functionT *) ParaM)->tableName="muscatt";
// fill table
((struct functionT *) ParaM)->mTable= gsl_matrix_alloc(nd, 19);
for ( int i=0; imTable,i,0,Qc[i] );
gsl_matrix_set( ((struct functionT *) ParaM)->mTable,i,1,III[i] );
gsl_matrix_set( ((struct functionT *) ParaM)->mTable,i,2,dIII[i] );
for ( int j=0; j<16; j++) gsl_matrix_set( ((struct functionT *) ParaM)->mTable,i,j+3,In[i][j]);
}
}
[code]
/* qtisas initial part: start */
int activeCol=int(activeOutput*2-2); // what to show after pressing "Simulate"
int numberPointsReal=currentLastPoint-currentFirstPoint+1;
if ( !beforeFit && !afterFit && !beforeIter && !afterIter) if (calculated)
{
if (activeCol>=0 && activeCol<15) return In[orderQ[currentPoint-currentFirstPoint]][activeCol];
else return YYY[currentPoint];
}
if ( afterIter) { calculated = false; return 1.0;};
if (currentInt!=0) { calculated = false; return 1.0;};
if (afterFit && calculated ) { saveToTable(numberPointsReal, ParaM ); return 1.0;};
/* qtisas initial part: finish */
if (numberPointsReal>300) numberPointsReal=300;
int orderQtmp[300];
double Qtab[300];
double Itab[300];
double Ieee[300];
int nm[300];
// double In[numberPoints][16]; =>> now defined as static
double AA[300];
int cnm;
double cQ,cI,cIe;
int i,j,k;
int nd;
double gR,gSl,gY0;
double pR,pSl,pY0;
double XX[100];
double YY[100];
double Ye[100];
double par[3];
double hhqcut;
double AAcorr,alpha;
double Qstep,Qmin;
double bgkerr,bgk,byy,Qcc;
int bnn,flag,steps;
double qv,its,wght;
//void ddct2d(int, int, int, double **, double *, int *, double *);
int *ip, n1, n2, n;
double **a, **b, *w;
double pi;
double thrd,sxth,svth;
int pos;
//
double I0;
double fac1,fac2,fac3,arg,eps;
double Tia,scl;
double sigmares,rr,resol,ddq;
double wt,twt,avg;
int num;
// para check
if (lambda<=0.0) {lambda=5.0;};
if (d<=0.0) {d=0.1;};
if (transm<=0.0 || transm>1.0) {transm=1.0;};
if (sigmia<=0.0) {sigmia=0.0;};
if (hgqmod<=0||hgqmod>1||hgqcut<=0.0) {hgqmod = 0;};
hhqcut = 0.0;
if (resmod<0 || resmod>6) {resmod = -1;}; // no resolution calculation
if (dlfwhm<0.0) {resmod = -1;};
if (dlfwhm>0.3) {dlfwhm = 0.1;}; // reset unreasonable high values
if (detlen<0.0) {resmod = -1;}
if (detlen>40.0) {detlen = 20.0;}; // reset unreasonable high values
if (sampap<=0.0) {resmod = -1;};
if (sampap>0.05) {sampap = 0.01;}; // reset unreasonable high values
if (collen<=0.0||entrap<=0.0){collen=-1.0; entrap=-1.0;}; // switch details about collimaion off
if (collen>40.0) {collen = detlen;};
if (entrap>0.10) {entrap = 2.0*sampap;};
if (reswdt<=0.0||reswdt>6.0) {reswdt = 2.0;};
if (thkmod<=0 || thkmod>12) {thkmod = -1;}; // no thickness correction
nd=numberPointsReal; //+++VP
for (i=0;i3||nd<3) lowQ = 3;
for (i=0;igR) {
gR = par[0]; // correlation coefficient
gSl= par[1]; // slope
gY0= par[2]; // ordinate value
};
};
};
// fprintf(stdout,"%12.4e %12.4e %12.4e \n",gR,gSl,gY0);
pR =-1.0;
pSl= 0.0;
pY0= 0.0;
if (lowQ==0||lowQ==2) {
for (i=0;ipR) {
pR = par[0]; // correlation coefficient
pSl= par[1]; // slope
pY0= par[2]; // ordinate value
};
};
};
// fprintf(stdout,"%12.4e %12.4e %12.4e \n",pR,pSl,pY0);
// only go on when more than 3 data points
if (nd>2) {
for (i=0;i<3;i++) {
XX[i] = (double)i;
YY[i] = Qc[i];
Ye[i] = 1.0;
};
linreg(XX,YY,Ye,i,par);
Qstep = par[1];
j = 3;
Ye[0] = 1.0;
Ye[1] = 1.0;
Ye[2] = 1.0;
YY[0] = 0.0; // slope before data begins
XX[0] = Qc[0]-0.5*Qstep;
if (lowQ<2) YY[0] = 2.0*gSl*XX[0]*exp(gY0+gSl*XX[0]*XX[0]);
if (lowQ==2||(lowQ=0 && pR>gR)) YY[0] = pSl*pow(XX[0],pSl-1.0)*exp(pY0);
YY[1] = (In[1][0]-In[0][0])/(Qc[1]-Qc[0]);
XX[1] = 0.5*(Qc[0]+Qc[1]);
for (i=0;i=0.87) {AA[i]=0.5*par[1];} else {AA[i]=0.0;};
YY[0] = YY[1];
YY[1] = YY[2];
XX[0] = XX[1];
XX[1] = XX[2];
};
bgkerr = 0.0;
if (bgksub==0.0) {
bgk = 0.0;
byy = 0.0;
bnn = 0;
for (i=nd-1;i>imax(0,nd-11);i--) {
bgk += In[i][0];
byy += In[i][0]*In[i][0];
bnn++;
};
bgksub = bgk / bnn;
bgkerr = sqrt(byy / bnn - bgksub*bgksub);
flag = 0;
while (flag==0) {
flag = 1;
if (i>=0) {
for (j=0;j=0) {
for (j=0;jQc[nd-1]) {hgqcut = Qc[i-j+1]*0.99999; hhqcut=Qc[i]*1.00001;};
Qcc = exp((log(0.5*bgkerr)-par[2])/par[1]);
// fprintf(stdout,"%12.4e\n", Qcc);
k = nd-1;
while (k>0 && Qc[k]>Qcc) {k--;};
if (nd-k<11) {k=imax((nd-11+k)/2,0);};
k++;
} else {k=0;};
// fprintf(stdout,"%i\n", k);
bgk = 0.0;
bnn = 0;
for (;khhqcut && i>0) {i--;};
i = imax(0,i-9);
hgqcut = Qc[i]*0.99999;
};
i = 0;
while(Qc[i]=3) {
linreg(XX,YY,Ye,j,par);
if (par[1]<-2.0) {
AAcorr = -0.5*pi/(2.0+par[1])*exp(par[2]);
alpha = par[1];
// fprintf(stdout,"%12.4e %12.4e\n",par[1],par[2]);
};
};
};
steps = (int)(pow(0.5,floor(log(Qc[nd-1]/Qstep)/log(0.5)))+0.5);
Qmin = 0.5*Qc[nd-1]/(steps-0.5);
Qstep = (Qc[nd-1]-Qmin)/(steps-1);
n1 = steps;
n2 = steps;
AAcorr *= (steps+0.5)*(steps+0.5)*pow(Qc[nd-1]+Qmin,alpha);
//fprintf(stdout,"%12.4e\n",AAcorr);
a = alloc_2d_double(n1, n2);
b = alloc_2d_double(n1, n2);
n = imax(n1, n2/2);
ip = alloc_1d_int(2 + (int)sqrt(n + 0.5));
n = imax(n1, n2) * 3 / 2;
w = alloc_1d_double(n);
// I0 = 0.0; total sum, i.e. field a(r=0)
for (i=0;igR)) its = pow(qv,pSl)*exp(pY0);
} else {
if (qv>Qc[nd-1]) {its=0.0;} else {
k=0;
while (Qc[k]=sampap/detlen) {
sigmares = (pow(Qc[k]*dlfwhm,2)
+ pow(2.0*pi/lambda*(entrap/collen-0.25*sampap*sampap/entrap*pow(detlen+collen,2)/(detlen*detlen*collen)),2) )
/ (8.0*log(2.0));
} else {
sigmares = (pow(Qc[k]*dlfwhm,2)
+ pow(2.0*pi/lambda*(entrap/collen+entrap/detlen-0.25*entrap*entrap/sampap*detlen/collen/(detlen+collen)),2) )
/ (8.0*log(2.0));
};
};
// fprintf(stdout,"%12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n",sigmares,entrap,collen,sampap,detlen,dlfwhm);
if (resmod==0) {
for (i=0;i0.0) {avg=avg/twt;};
In[k][pos+1]=avg*fac3;
};
};
// thickness correction
for (k=0;k6) {
if (collen<0.0) {
sigmares = (pow(Qc[k]*dlfwhm,2)
+ pow(2.0*pi/lambda*(1.5*sampap/detlen),2) ) / (8.0*log(2.0));
} else {
if (entrap/(collen+detlen)>=sampap/detlen) {
sigmares = (pow(Qc[k]*dlfwhm,2)
+ pow(2.0*pi/lambda*(entrap/collen-0.25*sampap*sampap/entrap*pow(detlen+collen,2)/(detlen*detlen*collen)),2) )
/ (8.0*log(2.0));
} else {
sigmares = (pow(Qc[k]*dlfwhm,2)
+ pow(2.0*pi/lambda*(entrap/collen+entrap/detlen-0.25*entrap*entrap/sampap*detlen/collen/(detlen+collen)),2) )
/ (8.0*log(2.0));
};
};
} else {
sigmares = 0.0;
};
// fprintf(stdout,"%12.4e %12.4e %12.4e %12.4e %i\n",fac1,fac2,sigmares,a[0][0],thkmod);
if ((thkmod>0 && thkmod<4) || (thkmod>6 && thkmod<10)) {
for (i=0;i6) {
for (i=0;i0.0) {avg=avg/twt;};
In[k][pos+1]=avg*fac3;
} else {
In[k][pos+1]=Ieee[nm[k]]*fac3;
};
};
};
/*
for (i=0;i=0 && activeCol<15) return In[orderQ[currentPoint-currentFirstPoint]][activeCol];
else return YYY[currentPoint];
//return In[orderQ[currentPoint-currentFirstPoint]][0];
[fortran]
0,
[end]