/******Selfing case******/
Data a;
  size=20;
  do itp=10 to 10;
    d1=itp*0.01; d2=(size-itp)*0.01; d=d1+d2;
     r1=0.5*(1-exp(-2.0*d1)) ; r2=0.5*(1-exp(-2.0*d2));
     r=0.5*(1-exp(-2.0*d));
     r1=0.1; r2=0.1; 
     p1=(1-r1)*(1-r2)/2;
     p2=r1*r2/2;
     p3=(1-r1)*r2/2;
     p4=r1*(1-r2)/2;
 
     A=p1**2; B=2*p1*p2; C=p2**2; D=2*p1*p3; E=2*p1*p4; F=2*p2*p3; G=2*p2*p4;
     H=p3**2; I=2*p3*p4; J=p4**2; K=2*p1*p4; L=2*p1*p3; M=2*p2*p4; N=2*p2*p3;
     S=2*p1*p2; T=2*p3*p4; U=2*p1**2; V=2*p2**2; W=2*p3**2; Z=2*p4**2;
     CC=(1-r)**2/4; DD=r**2/4; EE=r*(1-r)/2; FF=(1-r)**2/2; GG=r**2/2;
  
	 do iself=1 to 100;          
         
       A1=A+B/4+D/4+(1-r2)**2*E/4+r2**2*F/4+K/4+(1-r1)**2*L/4+r1**2*M/4
          +((1-r1)*(1-r2)+r1*r2)**2*S/4+(r1*(1-r2)+r2*(1-r1))**2*T/4
          +((1-r1)*(1-r2))**2*U/4+r1**2*r2**2*V/4+((1-r1)*r2)**2*W/4
          +(r1*(1-r2))**2*Z/4;
      B1=B/2+r2*(1-r2)*E/2+ r2*(1-r2)*F/2+r1*(1-r1)*L/2+ r1*(1-r1)*M/2
          +r1*r2*(1-r1)*(1-r2)*U/2 + r1*r2*(1-r1)*(1-r2)*V/2
          + r1*r2*(1-r1)*(1-r2)*W/2+ r1*r2*(1-r1)*(1-r2)*Z/2;
      C1=B/4+C+r2**2*E/4+(1-r2)**2*F/4+G/4+r1**2*L/4+(1-r1)**2*M/4
         +N/4+((1-r1)*(1-r2)+ r1*r2)**2*S/4+(r1*(1-r1)+r2*(1-r2))**2*T/4
         +( r1*r2)**2*U/4+((1-r1)*(1-r2))**2*V/4+(r1*(1-r2))**2*W/4
         +((1-r1)*r2)**2*Z/4;
      D1=D/2+r2*(1-r2)*E/2+r2*(1-r2)*F/2
         +((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2)+(1-r1)*r2)*S/2
         +((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2)+ (1-r1)*r2)*T/2
         +(1-r1)**2*r2*(1-r2)*U/2+r1**2*r2*(1-r2)*V/2
         +(1-r1)**2* r2*(1-r2)*W/2+r1**2* r2*(1-r2)*Z/2;
      E1=(1-r2)**2*E/2+r2**2*F/2+r1*(1-r1)*(1-r2)**2*U/2
          +r1*(1-r1)*r2**2*V/2+r1*(1-r1)*r2**2*W/2
          +r1*(1-r1)*(1-r2)**2*Z/2;
      F1=r2**2*E/2+(1-r2)**2*F/2+ r1*(1-r1)*r2**2*U/2+r1*(1-r1)*(1-r2)**2*V/2
         +r1*(1-r1)*(1-r2)**2*W/2+r1*(1-r1)*r2**2*Z/2;
      G1=r2*(1-r2)*E/2+r2*(1-r2)*F/2+G/2+((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2) +(1-r1)*r2)*S/2
         +((1-r1)*(1-r2)+ r1*r2)*(r1*(1-r2)+(1-r1)*r2)*T/2
         +r1**2*r2*(1-r2)*U/2+(1-r1)**2*r2*(1-r2)*V/2+r1**2* r2*(1-r2)*W/2
         +(1-r1)**2*r2*(1-r2)*Z/2;
      H1=D/4+r2**2*E/4+(1-r2)**2*F/4+H+(1-r1)**2*L/4+r1**2*M/4+I/4+N/4
         +(r1*(1-r2)+(1-r1)*r2)**2*S/4+((1-r1)*(1-r2)+r1*r2)**2*T/4
         +((1-r1)*r2)**2*U/4+(r1*(1-r2))**2*V/4+((1-r1)*(1-r2))**2*W/4
         +(r1*r2)**2*Z/4;
      I1=r2*(1-r2)*E/2+r2*(1-r2)*F/2+I/2+r1*(1-r1)*L/2+r1*(1-r1)*M/2
         +r1*(1-r1)*r2*(1-r2)*U/2+r1*(1-r1)*r2*(1-r2)*V/2
         +r1*(1-r1)*r2*(1-r2)*W/2+r1*(1-r1)*r2*(1-r2)*Z/2;
      J1=(1-r2)**2*E/4+r2**2*F/4+G/4+I/4+J+K/4+r1**2*L/4+(1-r1)**2*M/4
         +(r1*(1-r2)+(1-r1)*r2)**2*S/4+((1-r1)*(1-r2)+r1*r2)**2*T/4
         +(r1*(1-r2))**2*U/4+((1-r1)*r2)**2*V/4+(r1*r2)**2*W/4
         +((1-r1)*(1-r2))**2*Z/4;
      K1=K/2+r1*(1-r1)*L/2+r1*(1-r1)*M/2
          +((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2)+(1-r1)*r2)*S/2
          +((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2)+ (1-r1)*r2)*T/2
          +r1*(1-r1)*(1-r2)**2*U/2+r1*(1-r1)*r2**2*V/2
          +r1*(1-r1)*r2**2*W/2+r1*(1-r1)*(1-r2)**2*Z/2;
      L1=(1-r1)**2*L/2+r1**2*M/2+(1-r1)**2*r2*(1-r2)*U/2+r1**2*r2*(1-r2)*V/2
         +(1-r1)**2*r2*(1-r2)*W/2+r1**2*r2*(1-r2)*Z/2;
      M1=r1**2*L/2+(1-r1)**2*M/2+ r1**2*r2*(1-r2)*U/2+(1-r1)**2*r2*(1-r2)*V/2
         +r1**2*r2*(1-r2)*W/2+(1-r1)**2*r2*(1-r2)*Z/2;
      N1=r1*(1-r1)*L/2+r1*(1-r1)*M/2+N/2
         +((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2)+(1-r1)*r2)*S/2
         +((1-r1)*(1-r2)+r1*r2)*(r1*(1-r2)+(1-r1)*r2)*T/2
         +r1*(1-r1)*r2**2*U/2+r1*(1-r1)*(1-r2)**2*V/2
         +r1*(1-r1)*(1-r2)**2*W/2+r1*(1-r1)*r2**2*Z/2;
      S1=((1-r1)*(1-r2)+r1*r2)**2*S/2+(r1*(1-r2)+(1-r1)*r2)**2*T/2
          +r1*(1-r1)*r2*(1-r2)*U/2+r1*(1-r1)*r2*(1-r2)*V/2
          +r1*(1-r1)*r2*(1-r2)*W/2+r1*(1-r1)*r2*(1-r2)*Z/2;
      T1=(r1*(1-r2)+(1-r1)*r2)**2*S/2+((1-r1)*(1-r2)+ r1*r2)**2*T/2
          +r1*(1-r1)*r2*(1-r2)*U/2+r1*(1-r1)*r2*(1-r2)*V/2
          +r1*(1-r1)*r2*(1-r2)*W/2+r1*(1-r1)*r2*(1-r2)*Z/2;
      U1=((1-r1)*(1-r2))**2*U/2+(r1*r2)**2*V/2+((1-r1)*r2)**2*W/2
          +(r1*(1-r2))**2*Z/2;
      V1=(r1*r2)**2*U/2+((1-r1)*(1-r2))**2*V/2+(r1*(1-r2))**2*W/2
         +((1-r1)*r2)**2*Z/2;
      W1=((1-r1)*r2)**2*U/2+(r1*(1-r2))**2*V/2+((1-r1)*(1-r2))**2*W/2
         +(r1*r2)**2*Z/2;
      Z1=(r1*(1-r2))**2*U/2+((1-r1)*r2)**2*V/2+(r1*r2)**2*W/2
         +((1-r1)*(1-r2))**2*Z/2;
      sum=2*(A1+B1+C1+D1+E1+F1+G1+H1+I1+J1+K1+L1+M1+N1+S1+T1)+U1+V1+W1+Z1;

          A=A1; B=B1; C=C1; D=D1; E=E1; F=F1; G=G1;
          H=H1; I=I1; J=J1; L=L1; K=K1; M=M1; N=N1; S=S1;
          T=T1; U=U1; V=V1; W=W1; Z=Z1;
    
		 output;
end;
end;
keep   A1 B1 C1 D1 E1 F1 G1 H1 I1 J1 L1 K1 M1 N1 S1 T1 U1 V1 W1 Z1 sum;
run;
proc print data=a; run;




/******Random mating case******/

Data a;
  size=10;
  do itp=5 to 5;
  d1=itp*0.01; d2=(size-itp)*0.01; d=d1+d2;
   r1=0.5*(1-exp(-2.0*d1)) ; r2=0.5*(1-exp(-2.0*d2));    
   r=0.5*(1-exp(-2.0*d));
   r1=0.1; r2=0.1;
    p1=(1-r1)*(1-r2)/2;
   p2=r1*r2/2;
   p3=(1-r1)*r2/2;
   p4=r1*(1-r2)/2;
   p5=p4;
   p6=p3;
   p7=p2;
   p8=p1;
 do irandom=1 to 200;
	 
p1t=p1+r1*(p2*p4+p4*p5+p3*p4-p1*p6-p1*p7-p1*p8)+r2*(p2*p3+p3*p6+p3*p4-p1*p5-p1*p7-p1*p8)+r1*r2*(2*p1*p7+p1*p8+p2*p7-p3*p6-p4*p5-2*p3*p4);
p2t=p2+r1*(p1*p6+p3*p6+p5*p6-p2*p4-p2*p7-p2*p8)+r2*(p1*p5+p4*p5+p5*p6-p2*p3-p2*p7-p2*p8)+r1*r2*(p1*p8+p2*p7-p3*p6-p4*p5+2*p2*p8-2*p5*p6);
p3t=p3+r1*(p1*p7+p2*p7+p5*p7-p3*p4-p3*p6-p3*p8)+r2*(p1*p5+p1*p7+p1*p8-p2*p3-p3*p4-p3*p6)+r1*r2*(p3*p6+p4*p5-p1*p8-p2*p7+2*p3*p4-2*p1*p7);
p4t=p4+r1*(p1*p6+p1*p7+p1*p8-p2*p4-p3*p4-p4*p5)+r2*(p1*p7+p2*p7+p6*p7-p3*p4-p4*p5-p4*p8)+r1*r2*(p4*p5+p3*p6-p1*p8-p2*p7+2*p3*p4-2*p1*p7);
p5t=p5+r1*(p1*p8+p2*p8+p3*p8-p4*p5-p5*p6-p5*p7)+r2*(p2*p3+p2*p7+p2*p8-p1*p5-p4*p5-p5*p6)+r1*r2*(p3*p6+p4*p5-p1*p8-p2*p7+2*p5*p6-2*p2*p8);
p6t=p6+r1*(p2*p4+p2*p7+p2*p8-p1*p6-p3*p6-p5*p6)+r2*(p1*p8+p2*p8+p4*p8-p3*p6-p5*p6-p6*p7)+r1*r2*(p3*p6+p4*p5-p1*p8-p2*p7+2*p5*p6-2*p2*p8);
p7t=p7+r1*(p3*p4+p3*p6+p3*p8-p1*p7-p2*p7-p5*p7)+r2*(p3*p4+p4*p5+p4*p8-p1*p7-p2*p7-p6*p7)+r1*r2*(p1*p8+p2*p7-p3*p6-p4*p5+2*p1*p7-2*p3*p4);
p8t=p8+r1*(p4*p5+p5*p6+p5*p7-p1*p8-p2*p8-p3*p8)+r2*(p3*p6+p5*p6+p6*p7-p1*p8-p2*p8-p4*p8)+r1*r2*(p1*p8+p2*p7-p3*p6-p4*p5+2*p2*p8-2*p5*p6);

 p1=p1t;  p2=p2t;  p3=p3t;  p4=p4t;    p5=p5t;  p6=p6t;  p7=p7t;  p8=p8t;
 sum=p1t+p2t+p3t+p4t+p5t+p6t+p7t+p8t;  
 output;
 end;
end;
keep irandom itp  p1 p2 p3 p4 p5 p6 p7 p8 sum;
run;
proc print data=a; run;