public static double[][] calculateDensity() {
		
		int numGridLines = getNumGridlines();
		double[][] density = new double[numGridLines][numGridLines];
		
		double i1 = Settings.getMinimumRange();
		// Calculate un-normalized density for each point in domain, looping across X1
		for (int x1Ind = 0; x1Ind < numGridLines; x1Ind++)
			{
			i1 += Settings.discretization;
			int[] i1RelevantIndices = findRelevantKIndices(i1, Settings.startLevel);
			int k1Max = i1RelevantIndices[1];
			int k1Min = i1RelevantIndices[0];					
			
			// Cycle through relevant X1 translates for the line
			for (int k1Ind = k1Min; k1Ind < k1Max; k1Ind++) {
				
				double phi1Here = scalePhisHere.get(x1Ind).get(k1Ind - k1Min);
				
				// Loop across X2 dimension
				double i2 = Settings.getMinimumRange();
				for (int x2Ind = 0; x2Ind < numGridLines; x2Ind++)
				{
					i2 += Settings.discretization;
					int[] i2RelevantIndices = findRelevantKIndices(i2, Settings.startLevel);
					int k2Max = i2RelevantIndices[1];
					int k2Min = i2RelevantIndices[0];
										
					// Cycle through relevant X2 translates
					for (int k2Ind = k2Min; k2Ind < k2Max; k2Ind++) {
						
						double phi2Here = scalePhisHere.get(x2Ind).get(k2Ind - k2Min);
						
						density[x1Ind][x2Ind] += Transform.scalingCoefficients[k1Ind][k2Ind] *
								phi1Here*phi2Here;
					} // End cycling relevant X2 translates
				} // End cycling across X2 gridlines
			} // End cycling relevant X1 translates
		} // End cycling across X1 gridlines
		
		if (Settings.waveletFlag) {
			density = addWaveDensity(density);
		}
		
		//Normalize density
		density = normalizeDensity(density);
		
		return(density);
		
	}