source: sasmodels/explore/beta/testPolyDiseperseSphricalBeta.m @ db1d9d5

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since db1d9d5 was cdd676e, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

Update Yun's code to produce <F2> and P columns

  • Property mode set to 100644
File size: 1.7 KB
Line 
1close all
2clear all
3
4R0=20;       %mean radius. unit: Angstroms
5PDI=0.1;       %unit: relative polydispersity
6volF=0.1;
7contrast=6e-6;    %contrast,angstrom^-1, scattering length density difference
8
9Filename='testPolydisperseGaussianSphere.dat';
10
11Q=[0.001:0.001:0.8];
12
13N=10;
14
15deltaR=R0*PDI;
16SplitN=100;
17R=[R0/N:deltaR/SplitN:R0+5*deltaR];
18
19FGaussian=1/sqrt(2*pi*deltaR^2)*exp(-(R-R0).^2/2/deltaR^2);
20
21Ni=FGaussian*deltaR/SplitN;
22%trapz(R,FGaussian)
23%sum(FGaussian*deltaR/SplitN)
24sum(Ni)
25
26
27Vi=4*pi*R.^3/3;
28AverageVolume=sum(Ni.*Vi);
29
30
31IQN=zeros(size(Q));
32%IQ=sphereFormFactor(Q,R0);
33
34for i=1:length(R)
35    IQN=IQN+Ni(i)*Vi(i)^2*contrast^2*sphereFormFactor(Q,R(i));
36    %IQ=1*AverageVolume^2*contrast^2*sphereFormFactor(Q,R0);
37    %IQ=sphereFormFactor(Q,R(500));
38end
39
40IQ=volF/AverageVolume*IQN;
41
42IQ=IQ*1e8;  %convert from A^-1 to cm^-1
43loglog(Q,IQ,'b.-')
44
45FQ=zeros(size(Q));
46
47for i=1:length(R)
48    FQ=FQ+Ni(i)*Vi(i)*contrast*3*(sin(Q*R(i)) - Q*R(i).*cos(Q*R(i))) ./ (Q*R(i)).^3;
49end
50
51beta=FQ.^2./IQN;
52figure
53semilogx(Q,beta,'g*-');
54
55%stop
56
57NormIQN=IQN/IQN(1);
58
59FQ2=beta.*IQ;
60
61normQ=Q*R0*2;
62Sq=HSS_SQ(volF,normQ);
63
64figure
65loglog(Q,IQ,'bo-');
66figure
67plot(Q,beta,'g*-');
68figure
69plot(Q,Sq,'b.-');
70figure
71Sq_eff=1+beta.*(Sq-1);
72plot(Q,Sq_eff,'gd-');
73
74fileID=fopen(Filename,'w');
75
76FileString=strcat('Filename:',Filename, '**R0=', num2str(R0),'**PDI=', num2str(PDI),...
77    '**contrast',num2str(contrast), '**volF=',num2str(volF));
78Outputformat='File format : Q, <FQ>, <FQ^2>, PQ, betaQ, SQ, SQ_Eff';
79fprintf(fileID,'%s\n', FileString);
80fprintf(fileID,'%s\n', Outputformat);
81for i=1:length(Q)
82    %file format : Q, <FQ>, <FQ^2>, PQ, betaQ, SQ, SQ_Eff
83    fprintf(fileID,'%e\t%e\t%e\t%e\t%e\t%e\t%e\n',...
84        Q(i),FQ(i),IQN(i),IQ(i),beta(i),Sq(i),Sq_eff(i));
85end
86
87fclose(fileID)
88
Note: See TracBrowser for help on using the repository browser.