001 function [e,m,sigma2]= snplot(S,N,nr)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042 [n,m]=size(S);
043 if m>n, S=S'; end
044 [n,m]=size(S);
045 if m~=1, disp(' First argument must be nx1, program will terminate.'), end
046
047 [n,m]=size(N);
048 if m>n, N=N'; end
049 [n,m]=size(N);
050 if m~=1, disp(' Second argument must be nx1, program will terminate.'), end
051
052 A=[ones(length(S),1) log(S)];
053 beta=inv(A'*A)*A'*log(N); m=-beta(2); e=exp(-beta(1));
054 Q0=log(N)'*log(N)-beta'*A'*log(N);
055 sigma2=Q0/(length(log(S))-2);
056 k=exp(sigma2/2);
057 borderS=[min(S) max(S)];
058 borderN=[min(N) max(N)];
059
060 number_of_s=99;
061
062 if nargin==3
063 if nr<10
064 if nr==1 | nr==3
065 plot(S,N,'.','markersize',12)
066 title('SN-data')
067 elseif nr==2 | nr==4
068 s=borderS(1):(borderS(2)-borderS(1))/number_of_s:borderS(2);
069 n=k/e*s.^(-m);
070 plot(s,n,'-',S,N,'.','markersize',12)
071 axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
072 title('SN-data with estimated N(s)'),xlabel('s'), ylabel('N(s)')
073 end
074
075 axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
076 xlabel('s'), ylabel('N')
077
078 else
079 if nr==11 | nr==13
080 plot(N,S,'.','markersize',12)
081 title('SN-data')
082 elseif nr==12 | nr==14
083 s=borderS(1):(borderS(2)-borderS(1))/number_of_s:borderS(2);
084 n=k/e*s.^(-m);
085 plot(n,s,'-',N,S,'.','markersize',12)
086 axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
087 title('SN-data with estimated N(s)')
088 end
089
090 axis([0 1.1*borderN(2) [0.9 1.1].*borderS ]);
091 xlabel('N'), ylabel('s')
092 end
093
094 if mod(nr,10) > 2
095 h=gca;
096 set(h,'YScale','log','XScale','log');
097 end
098
099 end
100
101 return
102
103
104
105
106 if nargin==3
107 if nr==1
108 plot(S,N,'.','markersize',12)
109 axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
110 title('SN-data'),xlabel('s'), ylabel('N')
111 elseif nr==2
112 s=borderS(1):(borderS(2)-borderS(1))/number_of_s:borderS(2);
113 n=k/e*s.^(-m);
114 plot(s,n,'-',S,N,'.','markersize',12)
115 axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
116 title('SN-data with estimated N(s)'),xlabel('s'), ylabel('N(s)')
117 elseif nr==3
118 axis([[0.9 1.1].*log(borderS) 0 1.1*log(borderN(2))]);
119 plot(log(S),log(N),'.','markersize',12)
120 title('log(S)/log(N)-data'),xlabel('log(s)'), ylabel('log(N(s))')
121 elseif nr==4
122 y=beta(1)+beta(2)*log(borderS);
123 axis([[0.9 1.1].*log(borderS) 0 1.1*log(borderN(2))]);
124 plot(log(borderS),y,'-',log(S),log(N),'.','markersize',12)
125 title('log(S)/log(N)-data with estimated log(N(s))')
126 xlabel('log(s)'), ylabel('log(N(s))')
127 end
128 end
129
130