PRO multifrac_synth, im

sz_x=512
sz_y=512

im = fltarr(sz_x,sz_y)
im(*,*)=1
tvim,im


;calculate the probability density function
;use 10000 points as DMZ 'get_rid()' produces a number
;between 0 and 10000. This dictates the 0.01 seperation
prob_dist_den = fltarr(10000)
j=0L

FOR i=0.01,0.5,0.01 DO BEGIN & $
    prob_dist_den(j:(j-1+4*i*100)) = i & $
    j=max(where(prob_dist_den NE 0.))+1 & $
ENDFOR

FOR i=0.51,0.99,0.01 DO BEGIN & $
    prob_dist_den(j:(j-1+4*(1-i)*100)) = i & $
    j=max(where(prob_dist_den NE 0.))+1 & $
ENDFOR

n=0
WHILE ((2^(n)) LE 512) DO n=n+1
power=findgen(n)
rev_power=reverse(power)
   
FOR i=0,n-2 DO BEGIN
    print,i
    box_size=2.^(rev_power(i))
    num_boxes=2.^(2*(power(i)))
    num_points=sz_x/box_size

    ;work out the coordinates
    calc_coords,size(im),box_size,x_coords,y_coords

    ;j has to a LONG variable as it can get very large
    j=lon64arr(1)
    
     FOR j=0.,num_points-1 DO BEGIN
    	FOR k=0.,num_points-1 DO BEGIN
    	    splitim  = split_image( im(x_coords(j):x_coords(j+1)-1, y_coords(k):y_coords(k+1)-1), prob_dist_den) 
    	    im(x_coords(j):x_coords(j+1)-1, y_coords(k):y_coords(k+1)-1) = splitim
    	END  ;k loop
    END   ;j loop

tvim,im
print,total(im)
ENDFOR
END