Changeset 7ffa8196 in sasview for sansmodels/src/c_models


Ignore:
Timestamp:
Jan 5, 2012 10:25:27 AM (13 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
2d1b700
Parents:
c637521
Message:

refactor RPA model

Location:
sansmodels/src/c_models
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_models/models.hh

    rc637521 r7ffa8196  
    2626 
    2727using namespace std; 
    28  
    29  
    30  
    31  
    32  
    33 class RPAModel{ 
    34 public: 
    35         // Model parameters 
    36         Parameter lcase_n; 
    37         Parameter ba; 
    38         Parameter bb; 
    39         Parameter bc; 
    40         Parameter bd; 
    41  
    42         Parameter Kab; 
    43         Parameter Kac; 
    44         Parameter Kad; 
    45         Parameter Kbc; 
    46         Parameter Kbd; 
    47         Parameter Kcd; 
    48  
    49         Parameter scale; 
    50         Parameter background; 
    51  
    52         Parameter Na; 
    53         Parameter Phia; 
    54         Parameter va; 
    55         Parameter La; 
    56  
    57         Parameter Nb; 
    58         Parameter Phib; 
    59         Parameter vb; 
    60         Parameter Lb; 
    61  
    62         Parameter Nc; 
    63         Parameter Phic; 
    64         Parameter vc; 
    65         Parameter Lc; 
    66  
    67         Parameter Nd; 
    68         Parameter Phid; 
    69         Parameter vd; 
    70         Parameter Ld; 
    71  
    72         // Constructor 
    73         RPAModel(); 
    74  
    75         // Operators to get I(Q) 
    76         double operator()(double q); 
    77         double operator()(double qx, double qy); 
    78         double calculate_ER(); 
    79         double evaluate_rphi(double q, double phi); 
    80 }; 
    81  
    8228 
    8329class ReflModel{ 

  • sansmodels/src/c_models/rpa.cpp

    r67424cd r7ffa8196  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
     25#include "rpa.h" 
     26 
    2627using namespace std; 
    2728 
    28 extern "C" { 
    29         #include "rpa.h" 
     29static double rpa_kernel(double dp[], double q) { 
     30  int lCASE = dp[0]; 
     31 
     32  double Na,Nb,Nc,Nd,Nab,Nac,Nad,Nbc,Nbd,Ncd; 
     33  double Phia,Phib,Phic,Phid,Phiab,Phiac,Phiad; 
     34  double Phibc,Phibd,Phicd; 
     35  double va,vb,vc,vd,vab,vac,vad,vbc,vbd,vcd; 
     36  double m; 
     37  double ba,bb,bc,bd; 
     38  double Xa,Xb,Xc,Xd; 
     39  double Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad; 
     40  double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd; 
     41  double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd; 
     42  double S0da,S0db,S0dc,Pdd,S0dd; 
     43  double Kaa,Kab,Kac,Kad,Kba,Kbb,Kbc,Kbd; 
     44  double Kca,Kcb,Kcc,Kcd,Kda,Kdb,Kdc,Kdd; 
     45  double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc; 
     46  double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33; 
     47  double Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33; 
     48  double ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33; 
     49  double N11,N12,N13,N21,N22,N23,N31,N32,N33; 
     50  double M11,M12,M13,M21,M22,M23,M31,M32,M33; 
     51  double S11,S12,S13,S14,S21,S22,S23,S24; 
     52  double S31,S32,S33,S34,S41,S42,S43,S44; 
     53  double La,Lb,Lc,Ld,Lad,Lbd,Lcd,Nav,Intg; 
     54  double scale, background; 
     55 
     56  int ii=13;  //dp[ii<13] = fittable, else fixed 
     57 
     58  Na=1000.0; 
     59  Nb=1000.0; 
     60  Nc=1000.0; 
     61  Nd=1000.0;  //DEGREE OF POLYMERIZATION 
     62  Phia=0.25; 
     63  Phib=0.25; 
     64  Phic=0.25; 
     65  Phid=0.25 ; //VOL FRACTION 
     66  Kab=-0.0004; 
     67  Kac=-0.0004; 
     68  Kad=-0.0004;  //CHI PARAM 
     69  Kbc=-0.0004; 
     70  Kbd=-0.0004; 
     71  Kcd=-0.0004; 
     72  La=1.0e-12; 
     73  Lb=1.0e-12; 
     74  Lc=1.0e-12; 
     75  Ld=0.0e-12;   //SCATT. LENGTH 
     76  va=100.0; 
     77  vb=100.0; 
     78  vc=100.0; 
     79  vd=100.0  ; //SPECIFIC VOLUME 
     80  ba=5.0; 
     81  bb=5.0; 
     82  bc=5.0; 
     83  bd=5.0;   //SEGMENT LENGTH 
     84 
     85  //lCASE was shifted to N-1 from the original code 
     86  if (lCASE <= 1){ 
     87    Phia=0.0000001; 
     88    Phib=0.0000001; 
     89    Phic=0.5; 
     90    Phid=0.5; 
     91    Nc=dp[ii+8]; 
     92    Phic=dp[ii+9]; 
     93    vc=dp[ii+10]; 
     94    Lc=dp[ii+11]; 
     95    bc=dp[3]; 
     96    Nd=dp[ii+12]; 
     97    Phid=dp[ii+13]; 
     98    vd=dp[ii+14]; 
     99    Ld=dp[ii+15]; 
     100    bd=dp[4]; 
     101    Kcd=dp[10]; 
     102    scale=dp[11]; 
     103    background=dp[12]; 
     104  } 
     105  else if ((lCASE > 1) && (lCASE <= 4)){ 
     106    Phia=0.0000001; 
     107    Phib=0.333333; 
     108    Phic=0.333333; 
     109    Phid=0.333333; 
     110    Nb=dp[ii+4]; 
     111    Phib=dp[ii+5]; 
     112    vb=dp[ii+6]; 
     113    Lb=dp[ii+7]; 
     114    bb=dp[2]; 
     115    Nc=dp[ii+8]; 
     116    Phic=dp[ii+9]; 
     117    vc=dp[ii+10]; 
     118    Lc=dp[ii+11]; 
     119    bc=dp[3]; 
     120    Nd=dp[ii+12]; 
     121    Phid=dp[ii+13]; 
     122    vd=dp[ii+14]; 
     123    Ld=dp[ii+15]; 
     124    bd=dp[4]; 
     125    Kbc=dp[8]; 
     126    Kbd=dp[9]; 
     127    Kcd=dp[10]; 
     128    scale=dp[11]; 
     129    background=dp[12]; 
     130  } 
     131  else { 
     132    Phia=0.25; 
     133    Phib=0.25; 
     134    Phic=0.25; 
     135    Phid=0.25; 
     136    Na=dp[ii+0]; 
     137    Phia=dp[ii+1]; 
     138    va=dp[ii+2]; 
     139    La=dp[ii+3]; 
     140    ba=dp[1]; 
     141    Nb=dp[ii+4]; 
     142    Phib=dp[ii+5]; 
     143    vb=dp[ii+6]; 
     144    Lb=dp[ii+7]; 
     145    bb=dp[2]; 
     146    Nc=dp[ii+8]; 
     147    Phic=dp[ii+9]; 
     148    vc=dp[ii+10]; 
     149    Lc=dp[ii+11]; 
     150    bc=dp[3]; 
     151    Nd=dp[ii+12]; 
     152    Phid=dp[ii+13]; 
     153    vd=dp[ii+14]; 
     154    Ld=dp[ii+15]; 
     155    bd=dp[4]; 
     156    Kab=dp[5]; 
     157    Kac=dp[6]; 
     158    Kad=dp[7]; 
     159    Kbc=dp[8]; 
     160    Kbd=dp[9]; 
     161    Kcd=dp[10]; 
     162    scale=dp[11]; 
     163    background=dp[12]; 
     164  } 
     165 
     166  Nab=pow((Na*Nb),(0.5)); 
     167  Nac=pow((Na*Nc),(0.5)); 
     168  Nad=pow((Na*Nd),(0.5)); 
     169  Nbc=pow((Nb*Nc),(0.5)); 
     170  Nbd=pow((Nb*Nd),(0.5)); 
     171  Ncd=pow((Nc*Nd),(0.5)); 
     172 
     173  vab=pow((va*vb),(0.5)); 
     174  vac=pow((va*vc),(0.5)); 
     175  vad=pow((va*vd),(0.5)); 
     176  vbc=pow((vb*vc),(0.5)); 
     177  vbd=pow((vb*vd),(0.5)); 
     178  vcd=pow((vc*vd),(0.5)); 
     179 
     180  Phiab=pow((Phia*Phib),(0.5)); 
     181  Phiac=pow((Phia*Phic),(0.5)); 
     182  Phiad=pow((Phia*Phid),(0.5)); 
     183  Phibc=pow((Phib*Phic),(0.5)); 
     184  Phibd=pow((Phib*Phid),(0.5)); 
     185  Phicd=pow((Phic*Phid),(0.5)); 
     186 
     187  Xa=q*q*ba*ba*Na/6.0; 
     188  Xb=q*q*bb*bb*Nb/6.0; 
     189  Xc=q*q*bc*bb*Nc/6.0; 
     190  Xd=q*q*bd*bb*Nd/6.0; 
     191 
     192  Paa=2.0*(exp(-Xa)-1.0+Xa)/pow(Xa,2.0); 
     193  S0aa=Na*Phia*va*Paa; 
     194  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
     195  S0ab=(Phiab*vab*Nab)*Pab; 
     196  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
     197  S0ac=(Phiac*vac*Nac)*Pac; 
     198  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
     199  S0ad=(Phiad*vad*Nad)*Pad; 
     200 
     201  S0ba=S0ab; 
     202  Pbb=2.0*(exp(-Xb)-1.0+Xb)/pow(Xb,2.0); 
     203  S0bb=Nb*Phib*vb*Pbb; 
     204  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
     205  S0bc=(Phibc*vbc*Nbc)*Pbc; 
     206  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
     207  S0bd=(Phibd*vbd*Nbd)*Pbd; 
     208 
     209  S0ca=S0ac; 
     210  S0cb=S0bc; 
     211  Pcc=2.0*(exp(-Xc)-1.0+Xc)/pow(Xc,2.0); 
     212  S0cc=Nc*Phic*vc*Pcc; 
     213  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
     214  S0cd=(Phicd*vcd*Ncd)*Pcd; 
     215 
     216  S0da=S0ad; 
     217  S0db=S0bd; 
     218  S0dc=S0cd; 
     219  Pdd=2.0*(exp(-Xd)-1.0+Xd)/pow(Xd,2.0); 
     220  S0dd=Nd*Phid*vd*Pdd; 
     221  //lCASE was shifted to N-1 from the original code 
     222  switch(lCASE){ 
     223  case 0: 
     224    S0aa=0.000001; 
     225    S0ab=0.000002; 
     226    S0ac=0.000003; 
     227    S0ad=0.000004; 
     228    S0bb=0.000005; 
     229    S0bc=0.000006; 
     230    S0bd=0.000007; 
     231    S0cd=0.000008; 
     232    break; 
     233  case 1: 
     234    S0aa=0.000001; 
     235    S0ab=0.000002; 
     236    S0ac=0.000003; 
     237    S0ad=0.000004; 
     238    S0bb=0.000005; 
     239    S0bc=0.000006; 
     240    S0bd=0.000007; 
     241    break; 
     242  case 2: 
     243    S0aa=0.000001; 
     244    S0ab=0.000002; 
     245    S0ac=0.000003; 
     246    S0ad=0.000004; 
     247    S0bc=0.000005; 
     248    S0bd=0.000006; 
     249    S0cd=0.000007; 
     250    break; 
     251  case 3: 
     252    S0aa=0.000001; 
     253    S0ab=0.000002; 
     254    S0ac=0.000003; 
     255    S0ad=0.000004; 
     256    S0bc=0.000005; 
     257    S0bd=0.000006; 
     258    break; 
     259  case 4: 
     260    S0aa=0.000001; 
     261    S0ab=0.000002; 
     262    S0ac=0.000003; 
     263    S0ad=0.000004; 
     264    break; 
     265  case 5: 
     266    S0ab=0.000001; 
     267    S0ac=0.000002; 
     268    S0ad=0.000003; 
     269    S0bc=0.000004; 
     270    S0bd=0.000005; 
     271    S0cd=0.000006; 
     272    break; 
     273  case 6: 
     274    S0ab=0.000001; 
     275    S0ac=0.000002; 
     276    S0ad=0.000003; 
     277    S0bc=0.000004; 
     278    S0bd=0.000005; 
     279    break; 
     280  case 7: 
     281    S0ab=0.000001; 
     282    S0ac=0.000002; 
     283    S0ad=0.000003; 
     284    break; 
     285  case 8: 
     286    S0ac=0.000001; 
     287    S0ad=0.000002; 
     288    S0bc=0.000003; 
     289    S0bd=0.000004; 
     290    break; 
     291  default : //case 9: 
     292    break; 
     293  } 
     294  S0ba=S0ab; 
     295  S0ca=S0ac; 
     296  S0cb=S0bc; 
     297  S0da=S0ad; 
     298  S0db=S0bd; 
     299  S0dc=S0cd; 
     300 
     301  Kaa=0.0; 
     302  Kbb=0.0; 
     303  Kcc=0.0; 
     304  Kdd=0.0; 
     305 
     306  Kba=Kab; 
     307  Kca=Kac; 
     308  Kcb=Kbc; 
     309  Kda=Kad; 
     310  Kdb=Kbd; 
     311  Kdc=Kcd; 
     312 
     313  Zaa=Kaa-Kad-Kad; 
     314  Zab=Kab-Kad-Kbd; 
     315  Zac=Kac-Kad-Kcd; 
     316  Zba=Kba-Kbd-Kad; 
     317  Zbb=Kbb-Kbd-Kbd; 
     318  Zbc=Kbc-Kbd-Kcd; 
     319  Zca=Kca-Kcd-Kad; 
     320  Zcb=Kcb-Kcd-Kbd; 
     321  Zcc=Kcc-Kcd-Kcd; 
     322 
     323  DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc); 
     324 
     325  T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT; 
     326  T12= (S0ac*S0cb - S0ab*S0cc)/DenT; 
     327  T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT; 
     328  T21= (S0bc*S0ca - S0ba*S0cc)/DenT; 
     329  T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT; 
     330  T23= (S0ac*S0ba - S0aa*S0bc)/DenT; 
     331  T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT; 
     332  T32= (S0ab*S0ca - S0aa*S0cb)/DenT; 
     333  T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT; 
     334 
     335  Y1=T11*S0ad+T12*S0bd+T13*S0cd+1.0; 
     336  Y2=T21*S0ad+T22*S0bd+T23*S0cd+1.0; 
     337  Y3=T31*S0ad+T32*S0bd+T33*S0cd+1.0; 
     338 
     339  X11=Y1*Y1; 
     340  X12=Y1*Y2; 
     341  X13=Y1*Y3; 
     342  X21=Y2*Y1; 
     343  X22=Y2*Y2; 
     344  X23=Y2*Y3; 
     345  X31=Y3*Y1; 
     346  X32=Y3*Y2; 
     347  X33=Y3*Y3; 
     348 
     349  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
     350 
     351  m=1.0/(S0dd-ZZ); 
     352 
     353  N11=m*X11+Zaa; 
     354  N12=m*X12+Zab; 
     355  N13=m*X13+Zac; 
     356  N21=m*X21+Zba; 
     357  N22=m*X22+Zbb; 
     358  N23=m*X23+Zbc; 
     359  N31=m*X31+Zca; 
     360  N32=m*X32+Zcb; 
     361  N33=m*X33+Zcc; 
     362 
     363  M11= N11*S0aa + N12*S0ab + N13*S0ac; 
     364  M12= N11*S0ab + N12*S0bb + N13*S0bc; 
     365  M13= N11*S0ac + N12*S0bc + N13*S0cc; 
     366  M21= N21*S0aa + N22*S0ab + N23*S0ac; 
     367  M22= N21*S0ab + N22*S0bb + N23*S0bc; 
     368  M23= N21*S0ac + N22*S0bc + N23*S0cc; 
     369  M31= N31*S0aa + N32*S0ab + N33*S0ac; 
     370  M32= N31*S0ab + N32*S0bb + N33*S0bc; 
     371  M33= N31*S0ac + N32*S0bc + N33*S0cc; 
     372 
     373  DenQ1=1.0+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31; 
     374  DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33; 
     375  DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33; 
     376  DenQ=DenQ1+DenQ2+DenQ3; 
     377 
     378  Q11= (1.0 + M22-M23*M32 + M33 + M22*M33)/DenQ; 
     379  Q12= (-M12 + M13*M32 - M12*M33)/DenQ; 
     380  Q13= (-M13 - M13*M22 + M12*M23)/DenQ; 
     381  Q21= (-M21 + M23*M31 - M21*M33)/DenQ; 
     382  Q22= (1.0 + M11 - M13*M31 + M33 + M11*M33)/DenQ; 
     383  Q23= (M13*M21 - M23 - M11*M23)/DenQ; 
     384  Q31= (-M31 - M22*M31 + M21*M32)/DenQ; 
     385  Q32= (M12*M31 - M32 - M11*M32)/DenQ; 
     386  Q33= (1.0 + M11 - M12*M21 + M22 + M11*M22)/DenQ; 
     387 
     388  S11= Q11*S0aa + Q21*S0ab + Q31*S0ac; 
     389  S12= Q12*S0aa + Q22*S0ab + Q32*S0ac; 
     390  S13= Q13*S0aa + Q23*S0ab + Q33*S0ac; 
     391  S14=-S11-S12-S13; 
     392  S21= Q11*S0ba + Q21*S0bb + Q31*S0bc; 
     393  S22= Q12*S0ba + Q22*S0bb + Q32*S0bc; 
     394  S23= Q13*S0ba + Q23*S0bb + Q33*S0bc; 
     395  S24=-S21-S22-S23; 
     396  S31= Q11*S0ca + Q21*S0cb + Q31*S0cc; 
     397  S32= Q12*S0ca + Q22*S0cb + Q32*S0cc; 
     398  S33= Q13*S0ca + Q23*S0cb + Q33*S0cc; 
     399  S34=-S31-S32-S33; 
     400  S41=S14; 
     401  S42=S24; 
     402  S43=S34; 
     403  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
     404 
     405  Nav=6.022045e+23; 
     406  Lad=(La/va-Ld/vd)*sqrt(Nav); 
     407  Lbd=(Lb/vb-Ld/vd)*sqrt(Nav); 
     408  Lcd=(Lc/vc-Ld/vd)*sqrt(Nav); 
     409 
     410  Intg=pow(Lad,2.0)*S11+pow(Lbd,2.0)*S22+pow(Lcd,2.0)*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; 
     411 
     412  Intg *= scale; 
     413  Intg += background; 
     414 
     415  return Intg; 
     416 
    30417} 
    31418 
    32419RPAModel :: RPAModel() { 
    33         lcase_n = Parameter(0.0); 
    34         ba = Parameter(5.0); 
    35         bb = Parameter(5.0); 
    36         bc = Parameter(5.0); 
    37         bd = Parameter(5.0); 
    38  
    39     Kab = Parameter(-0.0004); 
    40     Kac = Parameter(-0.0004); 
    41     Kad = Parameter(-0.0004); 
    42     Kbc = Parameter(-0.0004); 
    43     Kbd = Parameter(-0.0004); 
    44     Kcd = Parameter(-0.0004); 
    45  
    46     scale = Parameter(1.0); 
    47     background = Parameter(0.0); 
    48  
    49     Na = Parameter(1000.0); 
    50     Phia = Parameter(0.25); 
    51     va = Parameter(100.0); 
    52     La = Parameter(1.0e-012); 
    53  
    54     Nb = Parameter(1000.0); 
    55     Phib = Parameter(0.25); 
    56     vb = Parameter(100.0); 
    57     Lb = Parameter(1.0e-012); 
    58  
    59     Nc = Parameter(1000.0); 
    60     Phic = Parameter(0.25); 
    61     vc = Parameter(100.0); 
    62     Lc = Parameter(1.0e-012); 
    63  
    64     Nd = Parameter(1000.0); 
    65     Phid = Parameter(0.25); 
    66     vd = Parameter(100.0); 
    67     Ld = Parameter(0.0e-012); 
     420  lcase_n = Parameter(0.0); 
     421  ba = Parameter(5.0); 
     422  bb = Parameter(5.0); 
     423  bc = Parameter(5.0); 
     424  bd = Parameter(5.0); 
     425 
     426  Kab = Parameter(-0.0004); 
     427  Kac = Parameter(-0.0004); 
     428  Kad = Parameter(-0.0004); 
     429  Kbc = Parameter(-0.0004); 
     430  Kbd = Parameter(-0.0004); 
     431  Kcd = Parameter(-0.0004); 
     432 
     433  scale = Parameter(1.0); 
     434  background = Parameter(0.0); 
     435 
     436  Na = Parameter(1000.0); 
     437  Phia = Parameter(0.25); 
     438  va = Parameter(100.0); 
     439  La = Parameter(1.0e-012); 
     440 
     441  Nb = Parameter(1000.0); 
     442  Phib = Parameter(0.25); 
     443  vb = Parameter(100.0); 
     444  Lb = Parameter(1.0e-012); 
     445 
     446  Nc = Parameter(1000.0); 
     447  Phic = Parameter(0.25); 
     448  vc = Parameter(100.0); 
     449  Lc = Parameter(1.0e-012); 
     450 
     451  Nd = Parameter(1000.0); 
     452  Phid = Parameter(0.25); 
     453  vd = Parameter(100.0); 
     454  Ld = Parameter(0.0e-012); 
    68455 
    69456} 
    … …  
    76463 */ 
    77464double RPAModel :: operator()(double q) { 
    78         double dp[29]; 
    79         // Fill parameter array for IGOR library 
    80         // Add the background after averaging 
    81         //Fittable parameters 
    82         dp[0] = lcase_n(); 
    83  
    84         dp[1] = ba(); 
    85         dp[2] = bb(); 
    86         dp[3] = bc(); 
    87         dp[4] = bd(); 
    88  
    89         dp[5] = Kab(); 
    90         dp[6] = Kac(); 
    91         dp[7] = Kad(); 
    92         dp[8] = Kbc(); 
    93         dp[9] = Kbd(); 
    94         dp[10] = Kcd(); 
    95  
    96         dp[11] = scale(); 
    97         dp[12] = background(); 
    98  
    99         //Fixed parameters 
    100         dp[13] = Na(); 
    101         dp[14] = Phia(); 
    102         dp[15] = va(); 
    103         dp[16] = La(); 
    104  
    105         dp[17] = Nb(); 
    106         dp[18] = Phib(); 
    107         dp[19] = vb(); 
    108         dp[20] = Lb(); 
    109  
    110         dp[21] = Nc(); 
    111         dp[22] = Phic(); 
    112         dp[23] = vc(); 
    113         dp[24] = Lc(); 
    114  
    115         dp[25] = Nd(); 
    116         dp[26] = Phid(); 
    117         dp[27] = vd(); 
    118         dp[28] = Ld(); 
    119  
    120         return rpa_kernel(dp,q); 
     465  double dp[29]; 
     466  // Fill parameter array for IGOR library 
     467  // Add the background after averaging 
     468  //Fittable parameters 
     469  dp[0] = lcase_n(); 
     470 
     471  dp[1] = ba(); 
     472  dp[2] = bb(); 
     473  dp[3] = bc(); 
     474  dp[4] = bd(); 
     475 
     476  dp[5] = Kab(); 
     477  dp[6] = Kac(); 
     478  dp[7] = Kad(); 
     479  dp[8] = Kbc(); 
     480  dp[9] = Kbd(); 
     481  dp[10] = Kcd(); 
     482 
     483  dp[11] = scale(); 
     484  dp[12] = background(); 
     485 
     486  //Fixed parameters 
     487  dp[13] = Na(); 
     488  dp[14] = Phia(); 
     489  dp[15] = va(); 
     490  dp[16] = La(); 
     491 
     492  dp[17] = Nb(); 
     493  dp[18] = Phib(); 
     494  dp[19] = vb(); 
     495  dp[20] = Lb(); 
     496 
     497  dp[21] = Nc(); 
     498  dp[22] = Phic(); 
     499  dp[23] = vc(); 
     500  dp[24] = Lc(); 
     501 
     502  dp[25] = Nd(); 
     503  dp[26] = Phid(); 
     504  dp[27] = vd(); 
     505  dp[28] = Ld(); 
     506 
     507  return rpa_kernel(dp,q); 
    121508} 
    122509 
    … …  
    128515 */ 
    129516double RPAModel :: operator()(double qx, double qy) { 
    130         double q = sqrt(qx*qx + qy*qy); 
    131         return (*this).operator()(q); 
     517  double q = sqrt(qx*qx + qy*qy); 
     518  return (*this).operator()(q); 
    132519} 
    133520 
    … …  
    140527 */ 
    141528double RPAModel :: evaluate_rphi(double q, double phi) { 
    142         return (*this).operator()(q); 
     529  return (*this).operator()(q); 
    143530} 
    144531 
    … …  
    148535 */ 
    149536double RPAModel :: calculate_ER() { 
    150 //NOT implemented!!! 
    151         return 0.0; 
     537  //NOT implemented!!! 
     538  return 0.0; 
    152539} 
Note: See TracChangeset for help on using the changeset viewer.