diff --git a/main.cpp b/main.cpp
index c73ddd6..be7b47d 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,4 +1,5 @@
 #include <iostream>
+#include <sstream>
 
 #include <GL/glut.h>
 
@@ -6,6 +7,7 @@
 #include <stim/parser/arguments.h>
 #include <stim/visualization/obj.h>
 #include <stim/visualization/gl_spharmonics.h>
+#include <stim/math/constants.h>
 
 #define theta_scale		0.01
 #define phi_scale		0.01
@@ -119,6 +121,60 @@ void mouse_drag(int x, int y){
 	glutPostRedisplay();
 }
 
+float uniformRandom()
+{
+	return (  (float)(rand()))/(  (float)(RAND_MAX)); 
+}
+
+std::vector<stim::vec3 <float> >
+sample_sphere(int num_samples, float radius = 1.0)
+{
+
+	float solidAngle = 2*stim::PI;	///Solid angle to sample over
+	float PHI[2], Z[2], range;	///Range of angles in cylinderical coordinates
+	PHI[0] = solidAngle/2;		///project the solid angle into spherical coords
+	PHI[1] = asin(0);		///
+	Z[0] = cos(PHI[0]);		///project the z into spherical coordinates
+	Z[1] = cos(PHI[1]);		///
+	range = Z[0] - Z[1];		///the range of all possible z values.
+
+	float z, theta, phi;		/// temporary individual
+
+	std::vector<stim::vec3<float> > samples;
+
+	//srand(time(NULL));			///set random seed
+	srand(100);			///set random seed
+
+	for(int i = 0; i < num_samples; i++)
+	{
+		z = uniformRandom()*range + Z[1];
+		theta = uniformRandom()*stim::TAU;
+		phi = acos(z);
+		stim::vec3<float> sph(1, theta, phi);
+		stim::vec3<float> cart = sph.sph2cart();
+		sph[0] *= radius;
+		samples.push_back(cart);
+	}
+	samples.push_back(stim::vec3<float>(0.,0.,1.));
+	samples.push_back(stim::vec3<float>(0.,1.0,0.));
+	samples.push_back(stim::vec3<float>(0.,-1.,0.));
+
+
+	std::stringstream name;
+      for(int i = 0; i < num_samples; i++)
+           name << samples[i].str() << std::endl;
+           name << samples[num_samples].str() << std::endl;
+           name << samples[num_samples+1].str() << std::endl;
+           name << samples[num_samples+2].str() << std::endl;
+	
+    
+      std::ofstream outFile;
+      outFile.open("New_Pos_Vectors.txt");
+      outFile << name.str().c_str();
+
+	return samples;
+}
+
 void process_arguments(int argc, char* argv[]){
 
 	args.add("help", "prints this help");
@@ -128,6 +184,7 @@ void process_arguments(int argc, char* argv[]){
 	args.add("obj", "approximates a geometric object given as a Wavefront OBJ file", "", "filename");
 	args.add("out", "filename for outputting spherical harmonics coefficients", "", "filename");
 	args.add("zaxis", "render the z-axis as a green line");
+	args.add("pdf", "outputs the PDF if an OBJ files is given");
 
 	//process the command line arguments
 	args.parse(argc, argv);
@@ -212,12 +269,17 @@ void process_arguments(int argc, char* argv[]){
 
 		std::string filename = args["obj"].as_string(0);
 		unsigned int l = args["obj"].as_int(1);
+		int p = args["obj"].as_int(2);
+		std::cout << p << std::endl;
+		std::vector<stim::vec3<float> > sphere = sample_sphere(p);
+		p = p+3;
 
 		//create an obj object
 		stim::obj<double> object(filename);
 
 		//get the centroid of the object
 		stim::vec<double> c = object.centroid();
+		c[0] = 0; c[1] = 0; c[2] = 0;
 
 		//get the number of vertices in the model
 		unsigned int nV = object.numV();
@@ -236,18 +298,59 @@ void process_arguments(int argc, char* argv[]){
 		//generate the spherical PDF
 		stim::spharmonics<double> P;
 		P.pdf(spherical, l, l);
+		std::vector<float> weights;		///array of weights
+		if(args["pdf"].is_set())
+		{
+//			S.pdf(spherical, l, l);
+			for(int i = 0; i < p; i++)	///for each point on the sphere.
+			{
+				float val = 0;		///value starts with 0
+				for(int j = 0; j < nV; j++)		///for each point on surface
+				{
+					stim::vec3<float> star(object.getV(j)[0] - c[0],
+						object.getV(j)[1] - c[1], 
+						object.getV(j)[2] - c[2]);		///center each point on the model 
+//					val += abs(star.dot(sphere[i])); 		///sum the dot product of the centered point and the sphere.
+					if(star.dot(sphere[i]) > 0)
+						val += pow(star.dot(sphere[i]),8); 		///sum the dot product of the centered point and the sphere.
+				}
+				weights.push_back(val);
+			}
+			
+			S.mcBegin(l,l);
+			for(int i = 0; i < p; i++)
+			{
+				if(sphere[i] == stim::vec3<float>(0., 0., 1.))
+				{
+					std::cout << i  << sphere[i] << " " << weights[i] << std::endl;
+				}
+				if(sphere[i] == stim::vec3<float>(0., 1., 0.))
+				{
+					std::cout << i << sphere[i] << " " << weights[i] << std::endl;
+				}
+				if(sphere[i] == stim::vec3<float>(0., -1., 0.))
+				{
+					std::cout << i <<  sphere[i] << " " << weights[i] << std::endl;
+				}
+				stim::vec3<float> sph = sphere[i].cart2sph();
+				S.mcSample(sph[1], sph[2], weights[i]);
+			}
+			S.mcEnd();
 
-		//begin Monte-Carlo sampling, using the model vertices as samples
-		S.mcBegin(l, l);
-		double theta, phi, fx, px;
-		for(unsigned int i = 0; i < nV; i++){
-			theta = spherical[i][1];
-			phi = spherical[i][2];
-			fx = spherical[i][0];
-			px = P(theta, phi);
-			S.mcSample(theta, phi, fx / px);
 		}
-		S.mcEnd();
+		else{ 
+			//begin Monte-Carlo sampling, using the model vertices as samples
+			S.mcBegin(l, l);
+			double theta, phi, fx, px;
+			for(unsigned int i = 0; i < nV; i++){
+				theta = spherical[i][1];
+				phi = spherical[i][2];
+				fx = spherical[i][0];
+				px = P(theta, phi);
+				S.mcSample(theta, phi, fx / px);
+			}
+			S.mcEnd();
+		}
 	}
 
 	//if the user specifies an SH basis function
--
libgit2 0.21.4