law.m

这是用于国土资源调查时

function law
%UNTITLED1 Summary of this function goes here
%  Detailed explanation goes here
code=input('Input the fucking codes');
[c,r]=size(code);
for i=1:c
    p(i)=code(i,1);%人口
    a(i)=code(i,2);%面积
    pm(i)=p(i)/a(i);%人口密度
end
%按人口密度从小到大排列
for i=1:c
    for j=i+1:c
        if pm(j)<pm(i)
            c1=pm(i);
            c2=p(i);
            c3=a(i);
            pm(i)=pm(j);
            p(i)=p(j);
            a(i)=a(j);
            pm(j)=c1;
            p(j)=c2;
            a(j)=c3;
        end
    end
end
popu=0;
area=0;
sigmac=0;
for i=1:c
    popu=popu+p(i);
    area=area+a(i);
end
for i=1:c
    la(i)=p(i)/popu;
    lb(i)=a(i)/area;
    sigmala(i+1)=0;
    sigmalb(i+1)=0;
    for j=1:1:i
    sigmala(i+1)=la(j)+sigmala(i+1);
    sigmalb(i+1)=lb(j)+sigmalb(i+1);
    end
    sigmac=sigmac+0.5*abs(lb(i)-la(i));
end  
sigmala(1)=0;
sigmalb(1)=0;
scatter(sigmalb,sigmala,'*')
ylabel('人口累计百分比')
xlabel('面积累计百分比')
title('劳伦兹曲线')
hold on
plot(sigmalb,sigmala,'-')
%xishu=polyfit(sigmalb,sigmala,2);
%sigmala1=polyval(xishu,sigmalb);
%plot(sigmalb,sigmala1,'-')
x=0:1;
y=x;
plot(x,y)
sigmaleft=0;
sigmaright=0;
for i=1:c-1
sigmaleft=sigmalb(i)*sigmala(i+1)+sigmaleft;
sigmaright=sigmalb(i+1)*sigmala(i)+sigmaright;
end
sigmaleft=sigmaleft+sigmalb(c)*sigmala(1);
sigmaright=sigmaright+sigmalb(1)*sigmala(c);
format short
C=sigmac
G=sigmaleft-sigmaright