Changeset ff85cc5 in sasview for sansmodels/src


Ignore:
Timestamp:
Aug 28, 2009 12:24:34 PM (15 years ago)
Author:
Jae Cho <jhjcho@…>
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:
a50ef5e
Parents:
e2afadf
Message:

revert to the original igor code

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sansmodels/src/libigor/libStructureFactor.c

    rd2e8e5d rff85cc5  
    1 /*      SimpleFit.c 
     1/*     SimpleFit.c 
    22 
    33A simplified project designed to act as a template for your curve fitting function. 
     
    77#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h 
    88#include "libStructureFactor.h" 
     9 
     10 
    911//Hard Sphere Structure Factor 
    1012// 
     
    1517        double calp,cbeta,cgam,prefac,c,vstruc; 
    1618        double r,phi,struc; 
    17          
     19 
    1820        r = dp[0]; 
    1921        phi = dp[1]; 
     
    2628        // 
    2729        //  calculate the structure factor 
    28         //      
     30        // 
    2931        a = 2.0*q*r; 
    3032        asq = a*a; 
     
    4042        vstruc = 1.0/(1.0-c); 
    4143        struc = vstruc; 
    42          
     44 
    4345        return(struc); 
    4446} 
     
    5456        double lam,lam2,test,mu,alpha,beta; 
    5557        double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 
    56          
     58 
    5759        qv= q; 
    5860        rad = dp[0]; 
     
    6062        eps = dp[2]; 
    6163        tau = dp[3]; 
    62          
     64 
    6365        eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps); 
    64          
     66 
    6567        sig = 2.0 * rad; 
    6668        aa = sig/(1.0 - eps); 
     
    114116        // 
    115117        sq = 1.0/(aq*aq +bq*bq); 
    116          
     118 
    117119        return(sq); 
    118120} 
     
    132134        double sigma,eta,eta2,eta3,eta4,etam1,etam14,alpha,beta,gamm; 
    133135        double x,sk,sk2,sk3,sk4,t1,t2,t3,t4,ck; 
    134          
     136 
    135137        x= q; 
    136          
     138 
    137139        req = dp[0]; 
    138140        phis = dp[1]; 
    139141        edibkb = dp[2]; 
    140142        lambda = dp[3]; 
    141          
     143 
    142144        sigma = req*2.; 
    143145        eta = phis; 
    144146        eta2 = eta*eta; 
    145147        eta3 = eta*eta2; 
    146         eta4 = eta*eta3;       
    147         etam1 = 1. - eta;  
     148        eta4 = eta*eta3; 
     149        etam1 = 1. - eta; 
    148150        etam14 = etam1*etam1*etam1*etam1; 
    149151        alpha = (  pow((1. + 2.*eta),2) + eta3*( eta-4.0 )  )/etam14; 
     
    152154        // 
    153155        //  calculate the structure factor 
    154          
     156 
    155157        sk = x*sigma; 
    156158        sk2 = sk*sk; 
     
    165167        ck =  -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3; 
    166168        struc  = 1./(1.-ck); 
    167          
     169 
    168170        return(struc); 
    169171} 
    170172 
    171 // Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions  
     173// Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions 
    172174// 
    173175double 
     
    182184        double pi; 
    183185        int ierr; 
    184          
     186 
    185187        pi = 4.0*atan(1.); 
    186188        QQ= q; 
    187          
    188         diam=dp[0]*2;           //in   dp[0] coming from python is radius : cahnged on Mar. 30, 2009 
     189 
     190        diam=dp[0];             //in � 
    189191        zz = dp[1];             //# of charges 
    190         VolFrac=dp[2];   
     192        VolFrac=dp[2]; 
    191193        Temp=dp[3];             //in degrees Kelvin 
    192194        csalt=dp[4];            //in molarity 
     
    202204        Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0);   //in m^3 
    203205        cs=csalt*6.022E23*1E3;  //# salt molecules/m^3 
    204          
     206 
    205207        //         Compute the derived values of : 
    206         //                       Ionic strength IonSt (in C^2/m^3)   
     208        //                       Ionic strength IonSt (in C^2/m^3) 
    207209        //                      Kappa (Debye-Huckel screening length in m) 
    208210        //      and             gamma Exp(-k) 
     
    211213                                                                           //   Kappa=2/SIdiam                                  // Use to compare with HP paper 
    212214        gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2)); 
    213          
    214         //         Finally set up dimensionless parameters  
     215 
     216        //         Finally set up dimensionless parameters 
    215217        Qdiam=QQ*diam; 
    216218        gMSAWave[6] = Kappa*SIdiam; 
    217219        gMSAWave[4] = VolFrac; 
    218          
     220 
    219221        //Function sqhpa(qq)  {this is where Hayter-Penfold program began} 
    220          
     222 
    221223        //       FIRST CALCULATE COUPLING 
    222          
     224 
    223225        ss=pow(gMSAWave[4],(1.0/3.0)); 
    224226        gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss); 
    225          
     227 
    226228        //        CALCULATE COEFFICIENTS, CHECK ALL IS WELL 
    227229        //        AND IF SO CALCULATE S(Q*SIG) 
    228          
     230 
    229231        ierr=0; 
    230232        ierr=sqcoef(ierr); 
     
    237239                //      print "Please report HPMSA problem with above error code" 
    238240        } 
    239          
     241 
    240242        return(SofQ); 
    241243} 
     
    267269int 
    268270sqcoef(int ir) 
    269 {        
     271{ 
    270272        int itm=40,ix,ig,ii; 
    271273        double acc=5.0E-6,del,e1,e2,f1,f2; 
    272          
     274 
    273275        //      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" 
    274276        f1=0;           //these were never properly initialized... 
    275277        f2=0; 
    276          
     278 
    277279        ig=1; 
    278280        if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4])) { 
     
    290292        gMSAWave[10]=fmin(gMSAWave[4],0.20); 
    291293        if ((ig!=1) || ( gMSAWave[9]>=0.15)) { 
    292                 ii=0;                              
     294                ii=0; 
    293295                do { 
    294296                        ii=ii+1; 
    295297                        if(ii>itm) { 
    296298                                ir=-1; 
    297                                 return ir;               
     299                                return ir; 
    298300                        } 
    299301                        if (gMSAWave[10]<=0.0) { 
     
    347349int 
    348350sqfun(int ix, int ir) 
    349 {        
     351{ 
    350352        double acc=1.0e-6; 
    351353        double reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,rgek; 
     
    365367        p2=0; 
    366368        p3=0; 
    367          
     369 
    368370        //     CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981) 
    369          
    370         reta = gMSAWave[16];                                                 
     371 
     372        reta = gMSAWave[16]; 
    371373        eta2 = reta*reta; 
    372374        eta3 = eta2*reta; 
     
    379381                ibig=1; 
    380382        } 
    381      
     383 
    382384        gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12]); 
    383385        rgek =  gMSAWave[11]; 
     
    389391        d = 1.0-reta; 
    390392        d2 = d*d; 
    391         dak = d/rak;                                                   
     393        dak = d/rak; 
    392394        dd2 = 1.0/d2; 
    393395        dd4 = dd2*dd2; 
     
    400402        eta21 = 2.0*reta+1.0; 
    401403        eta22 = eta21*eta21; 
    402          
     404 
    403405        //     ALPHA(I) 
    404          
     406 
    405407        al1 = -eta21*dak; 
    406408        al2 = (14.0*eta2-4.0*reta-1.0)*dak2; 
    407409        al3 = 36.0*eta2*dak4; 
    408          
     410 
    409411        //      BETA(I) 
    410          
     412 
    411413        be1 = -(eta2+7.0*reta+1.0)*dak; 
    412414        be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2; 
    413415        be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4; 
    414          
     416 
    415417        //      NU(I) 
    416          
     418 
    417419        vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak; 
    418420        vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2; 
     
    420422        vu4 = vu1+e24*rak*vu3; 
    421423        vu5 = eta6d*(vu2+4.0*vu3); 
    422          
     424 
    423425        //      PHI(I) 
    424          
     426 
    425427        ph1 = eta6d/rak; 
    426428        ph2 = d-e12*dak2; 
    427          
     429 
    428430        //      TAU(I) 
    429          
     431 
    430432        ta1 = (reta+5.0)/(5.0*rak); 
    431433        ta2 = eta2d*dak2; 
     
    433435        ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2); 
    434436        ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2; 
    435          
     437 
    436438        //     double PRECISION SINH(K), COSH(K) 
    437          
     439 
    438440        ex1 = exp(rak); 
    439441        ex2 = 0.0; 
     
    445447        ckma = ck-1.0-rak*sk; 
    446448        skma = sk-rak*ck; 
    447          
     449 
    448450        //      a(I) 
    449          
     451 
    450452        a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4; 
    451453        if (ibig==0) { 
     
    453455                a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4; 
    454456        } 
    455          
     457 
    456458        //      b(I) 
    457          
     459 
    458460        b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4; 
    459461        if (ibig==0) { 
     
    461463                b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4; 
    462464        } 
    463          
     465 
    464466        //      V(I) 
    465          
     467 
    466468        v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45; 
    467469        if (ibig==0) { 
     
    469471                v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45; 
    470472        } 
    471          
    472          
     473 
     474 
    473475        //       P(I) 
    474          
     476 
    475477        pp1 = ph1*ph1; 
    476478        pp2 = ph2*ph2; 
     
    482484                p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2; 
    483485        } 
    484          
     486 
    485487        //       T(I) 
    486          
     488 
    487489        t1 = ta3+ta4*a1+ta5*b1; 
    488490        if (ibig!=0) { 
    489                  
     491 
    490492                //              VERY LARGE SCREENING:  ASYMPTOTIC SOLUTION 
    491                  
     493 
    492494                v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45; 
    493495                t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0; 
     
    516518                } 
    517519                gMSAWave[10] = gMSAWave[16]; 
    518                  
     520 
    519521        } else { 
    520          
     522 
    521523                t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk); 
    522524                t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0; 
    523                  
     525 
    524526                //              MU(i) 
    525                  
     527 
    526528                um1 = t2*a2-e12*v2*v2; 
    527529                um2 = t1*a2+t2*a1-e24*v1*v2; 
     
    530532                um5 = t1*a3+t3*a1-e24*v1*v3; 
    531533                um6 = t3*a3-e12*v3*v3; 
    532                  
     534 
    533535                //                      GILLAN CONDITION ? 
    534536                // 
     
    538540                // 
    539541                if ((ix==1) || (ix==3)) { 
    540                          
     542 
    541543                        //                      NO - CALCULATE REMAINING COEFFICIENTS. 
    542                          
     544 
    543545                        //                      LAMBDA(I) 
    544                          
     546 
    545547                        al1 = e12*p2*p2; 
    546548                        al2 = e24*p1*p2-b2-b2; 
     
    549551                        al5 = e24*p1*p3-b3-b3-ak2; 
    550552                        al6 = e12*p3*p3; 
    551                          
     553 
    552554                        //                      OMEGA(I) 
    553                          
     555 
    554556                        w16 = um1*al6-al1*um6; 
    555557                        w15 = um1*al5-al1*um5; 
     
    557559                        w13 = um1*al3-al1*um3; 
    558560                        w12 = um1*al2-al1*um2; 
    559                          
     561 
    560562                        w26 = um2*al6-al2*um6; 
    561563                        w25 = um2*al5-al2*um5; 
    562564                        w24 = um2*al4-al2*um4; 
    563                          
     565 
    564566                        w36 = um3*al6-al3*um6; 
    565567                        w35 = um3*al5-al3*um5; 
     
    571573                        w3526 = w35+w26; 
    572574                        w3425 = w34+w25; 
    573                          
     575 
    574576                        //                      QUARTIC COEFFICIENTS 
    575                          
     577 
    576578                        w4 = w16*w16-w13*w36; 
    577579                        w3 = 2.0*w16*w15-w13*w3526-w12*w36; 
     
    579581                        w1 = 2.0*w15*w14-w13*w24-w12*w3425; 
    580582                        w0 = w14*w14-w12*w24; 
    581                          
     583 
    582584                        //                      ESTIMATE THE STARTING VALUE OF f 
    583                          
     585 
    584586                        if (ix==1) { 
    585587                                //                              LARGE K 
     
    600602                                fap = -(pg+p2*ca)/p3; 
    601603                        } 
    602                          
     604 
    603605                        //                      AND REFINE IT ACCORDING TO NEWTON 
    604606                        ii=0; 
     
    616618                                del=fabs((fap-fa)/fa); 
    617619                        } while (del>acc); 
    618                          
     620 
    619621                        ir = ir+ii; 
    620622                        fa = fap; 
     
    648650} 
    649651 
    650 // called as DiamCyl(hcyl,rcyl)  
    651 //modified from Igor NIST package XOP 
     652// called as DiamCyl(hcyl,rcyl) 
    652653double 
    653 DiamCyl(double dp[], double q) 
    654 { 
    655         double hcyl, rcyl; 
     654DiamCyl(double hcyl, double rcyl) 
     655{ 
     656 
    656657        double diam,a,b,t1,t2,ddd; 
    657658        double pi; 
    658          
    659         rcyl = dp[0]; 
    660         hcyl = dp[1]; 
     659 
    661660        pi = 4.0*atan(1.0); 
    662         if (rcyl == 0 || hcyl == 0) { 
    663                 return 0.0; 
    664         } 
     661 
    665662        a = rcyl; 
    666663        b = hcyl/2.0; 
     
    669666        ddd = 3.0*t1*t2; 
    670667        diam = pow(ddd,(1.0/3.0)); 
    671          
    672         return(diam/2);  //return radius 
     668 
     669        return(diam); 
    673670} 
    674671 
     
    680677//returns DIAMETER 
    681678// called as DiamEllip(aa,bb) 
    682  
    683 //modified from Igor NIST package XOP 
    684679double 
    685 DiamEllip(double dp[], double q) 
    686 { 
    687         double aa, bb; 
     680DiamEllip(double aa, double bb) 
     681{ 
     682 
    688683        double ee,e1,bd,b1,bL,b2,del,ddd,diam; 
    689          
    690         aa = dp[0]; 
    691         bb = dp[1]; 
    692         if (aa == 0 || bb == 0) { 
    693                 return 0.0; 
    694         } 
    695         if (aa == bb) { 
    696                 return aa; 
    697         } 
     684 
    698685        if(aa>bb) { 
    699686                ee = (aa*aa - bb*bb)/(aa*aa); 
     
    701688                ee = (bb*bb - aa*aa)/(bb*bb); 
    702689        } 
    703          
     690 
    704691        bd = 1.0-ee; 
    705692        e1 = sqrt(ee); 
     
    708695        b2 = 1.0 + bd/2/e1*log(bL); 
    709696        del = 0.75*b1*b2; 
    710          
     697 
    711698        ddd = 2.0*(del+1.0)*aa*bb*bb;           //volume is always calculated correctly 
    712699        diam = pow(ddd,(1.0/3.0)); 
    713          
    714         return(diam/2);      //return radius 
     700 
     701        return(diam); 
    715702} 
    716703 
    717704double 
    718705sqhcal(double qq) 
    719 {        
    720     double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter;            
     706{ 
     707    double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter; 
    721708        //      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" 
    722          
     709 
    723710        etaz = gMSAWave[10]; 
    724711        akz =  gMSAWave[12]; 
     
    733720        sk = 0.5*(x1-x2); 
    734721        ak2 = akz*akz; 
    735          
     722 
    736723        if (qq<=0.0) { 
    737724                SofQ = -1.0/gMSAWave[0]; 
Note: See TracChangeset for help on using the changeset viewer.