*! version 1.0.0   01june1995
program define scatt
        version 4.0        
        
        preserve
        global HJN_OV1
        global HJN_OV2

        * To assign a default value include the following line
        global HJN_V1 = 100
        global HJN_1 "Sample Size"
        wdctl static HJN_1   15 10  55  8
        wdctl edit   HJN_V1  72 10  25  8 

        global HJN_V2 = 1
        global HJN_2 "Symbol Number"
        wdctl static HJN_2   15 25  55  8
        wdctl edit   HJN_V2  72 25  25  8 

        
        wdctl button "Run"   15 90 25 12 HJN_B2
        wdctl button "Close"  45 90 25 12 HJN_B3
        wdctl button "Help"  75 90 25 12 HJN_B4
        
        global HJN_B2 "scattdr"
        global HJN_B4 "whelp scatt"
        global HJN_B3 "exit 1234"
        
        noi wdlg "Scatterplots"  5 5 130 132 
        restore
end

program define scattdr

        preserve
        clear

        tempvar e1 e2

        scalar n=$HJN_V1

        if n< 2 | n > 600 {
                sstopbox stop "sample size must be between 2 and 600"
                        exit
                }


        local sym=$HJN_V2
        local nn=n
        quietly set obs `nn'

        scalar i = 1
        scalar j = 1
        scalar nc = 1000
        scalar nr = 6000
        scalar np = 4000
        scalar nix = 5000
        scalar niy = 5000
        scalar ii=1

        scalar rho = -1.0
        local mrho=rho

        gph open
        gph pen 1

        gen `e1' = uniform()
        gen `e2' = uniform()

        while i < 5 {

                scalar j = 1
                scalar nc = 1000

                while j < 6 {

                        if ii != 11 { scalar rho = rho + 0.1 }

                        local y1=nr-np
                        local y2=nr
                        local x1=nc
                        local x2=nc+np

                        gph line `y1' `x1' `y1' `x2' 1
                        gph line `y1' `x2' `y2' `x2' 1
                        gph line `y2' `x2' `y2' `x1' 1
                        gph line `y2' `x1' `y1' `x1' 1
                        
*                        gph box `y1' `x1' `y2' `x2' 4

                        replace `e1'=invnorm(uniform())
                        replace `e2'=rho*`e1'+sqrt(1.0-rho*rho)*invnorm(uniform())

                        replace `e1' = nc + np*((`e1' + 3)/6)
                        replace `e2' = nr - np*((`e2' + 3)/6)
                        gph vpoint `e2' `e1' , size(10) symbol(`sym')

                        local mrho: display %4.2f rho
                        local my=int(nr-np)-300
                        local mx=int(nc + (np/2))
                        gph text `my' `mx' 0 0 rho=`mrho' 

                        scalar j = j + 1
                        scalar nc = nc + nix

                        scalar ii=ii+1

                }

        scalar i = i + 1
        scalar nr = nr + niy

}

        gph close

        restore

end

exit