#######################################
## R function to draw a pps sample   ##
## Rao-Sampford pps sampling method  ##
##                                   ##
## Input: p is the Nx1 popu vector   ##
##        for the size variable      ##
##        n is the sample size       ##
## Output: the set of sampled units  ##
##                                   ##
## Written by Changbao Wu, July 2004 ##
#######################################

RSsample<-function(p,n){
N<-length(p)
lam<-rep(0,N+1)
for(i in 1:N) lam[i+1]<-lam[i]+p[i]
q<-rep(0,N)
for(i in 1:N) q[i]<-p[i]/(1-n*p[i])
q<-q/sum(q)
lam2<-rep(0,N+1)
for(i in 1:N) lam2[i+1]<-lam2[i]+q[i]
ntot<-1
while(ntot<1963){
sam<-NULL
rand<-runif(1)
dif<-1
L<-1
U<-N+1
while(dif>0){
M<-floor((U-L)/2)
if(lam[L+M]>rand) U<-U-M
if(lam[L+M]<=rand) L<-L+M
if(lam[U]>=rand & lam[U-1]<rand){
si<-U-1
dif<-0
      }
if(lam[L]<rand & lam[L+1]>=rand){
si<-L
dif<-0
      }
     }
sam<-c(sam,si)
nn<-0
while(nn<n-1){
dif<-1
rand<-runif(1)
L<-1
U<-N+1
while(dif>0){
M<-floor((U-L)/2)
if(lam2[L+M]>rand) U<-U-M
if(lam2[L+M]<=rand) L<-L+M
if(lam2[U]>=rand & lam2[U-1]<rand){
si<-U-1
dif<-0
      }
if(lam2[L]<rand & lam2[L+1]>=rand){
si<-L
dif<-0
      }
     }
if(min(abs(si-sam))==0) nn<-n+1
if(min(abs(si-sam))>0){
sam<-c(sam,si)
nn<-nn+1
          }
         }
if(nn==n+1) ntot<-1
if(nn==n-1) ntot<-1989
             }
return(sam)
           }