// test Hammersley Point Set in Maya

proc float[] hammersleyPointSet( int $pointNumber )
{
   float $pointSet[] = {};

   // find the num of decimals (to be replaced by a log base 2)
   int $numOfDecimals = 0;
   float $powerOfTwo = 0;
   while ( $powerOfTwo < $pointNumber )
   {
       $numOfDecimals += 1;
       $powerOfTwo = pow(2,$numOfDecimals);
   }
   
??// ratio between num of point and the next highest power of 2
   float $powOfTwoRatio = $powerOfTwo /float( $pointNumber );
   int $iter = 0;
   for( $iter = 0; $iter < $pointNumber; $iter ++ )
   {
             // GENERATE BINARY FRACTIONS
       float $binaryValues[] = { 0 };
       float $revBinValues[] = { 0 };

       int $decimal = 0;
       float $remainder = $iter;
       for( $decimal = ($numOfDecimals - 1); $decimal >= 0; $decimal -- )
       {
           float $power = pow( 2, $decimal );
           if ( $remainder >= $power  )
           {
               $binaryValues[ ($numOfDecimals-1) - $decimal ] = 1;
               $revBinValues[ $decimal ] = 1;
               $remainder -= $power;
           }
           else
           {
               $binaryValues[ ($numOfDecimals-1) - $decimal ] = 0;
               $revBinValues[ $decimal ] = 0;
           }

       }                
????
	// CONVERT BACK TO FLOATS
       float $x = 0;
       float $y = 0;
       int $index = 0;
       for( $index = 1; $index <=  $numOfDecimals; $index ++ )
       {
           $x += ( pow( 2, -$index ) * $binaryValues[ $index-1 ] );
           $y += ( pow( 2, -$index ) * $revBinValues[ $index-1 ] );
       }
       // scale to compensate to the "non-power-of-two" number of points
       $x *= $powOfTwoRatio;
        
?????????????
       $pointSet[ ($iter * 2) ] = $x;
       $pointSet[ ($iter * 2) + 1 ] = $y;
   }
       return $pointSet;
}

proc testHPS( int $samples )
{
	float $sampleCoord[] = hammersleyPointSet( $samples );
	string $cmd = "particle ";
	int $i = 0;
	for ( $i = 0; $i < ($samples * 2);  $i +=2 )
	{
		float $angleRad = $sampleCoord[$i] * 2 * 3.14159; // normalised to 2PI
		float $x =  cos( $angleRad ) * $sampleCoord[$i+1]; 
		float $y = sqrt( pow( $sampleCoord[$i+1], 2) - pow( $x, 2 ) );
		if( $angleRad > 3.14159 ) $y *= -1;
		$cmd += ("-p "+$x+" 0 "+$y+" ");
	}
	eval( $cmd );
}

proc HamersleyTest( int $samples, int $uniformSampling )
{
	float $sampleCoord[] = hammersleyPointSet( $samples );
	string $cmd = "particle ";
	int $i = 0;
	for ( $i = 0; $i < ($samples * 2);  $i +=2 )
	{
		float $coordX = $sampleCoord[$i];
		float $coordY = $sampleCoord[$i+1];
		if ( $uniformSampling ) $coordY = sqrt( $coordY );
		float $angleRad = $coordX * 2 * 3.14159; // normalised to 2PI
		float $x =  cos( $angleRad ) * $coordY; 
		float $y = sqrt( pow( $coordY, 2) - pow( $x, 2 ) );
		if( $angleRad > 3.14159 ) $y *= -1;
		$cmd += ("-p "+$x+" 0 "+$y+" ");
	}
	eval( $cmd );
}

HamersleyTest( 128, 0 );