# Integration in 3D

So the question was to use a nested loop to solve a 3D integral with the function conditions (written below in the code) to find $$\int dxdydz$$ and the x and y coordinate of the centre of mass of this domain. I am using 10 integrand evaluations per direction, but I am getting a large error especially while calculating the y center of mass coordinate. The limits of x were from 0 to 4, y from -4 to 4, z from -1 to 1. (And the density is 1) This is the outline of the code I wrote,

#include <stdio.h>
#include <math.h>

double func(double x, double y, double z){
double out;
double condition;
condition = z*z + pow((sqrt(x*x + y*y) - 3), 2); //defining the domain

if (condition <= 1 && x >= 0 && y >= -4){
out = 1;
}
else{
out = 0;
}

return out;
}

int main(){

int i,j,k;

double sum1 = 0, sum2 = 0;

for (i = 0; i < 10; i++){
double x = 0 + i*(4)/10;
double sum_y1 =0, sum_y2 =0;

for (j=0;j<10;j++){
double y = -4 + j*(8)/10;
double sum_z1 =0, sum_z2 =0;

for (k=0;k<10;k++){
double z = 1 - k*(2)/10;
double x_cm_value = x*func(x,y,z);

sum_z1 =  sum_z1 + func(x,y,z)*2/10;
sum_z2 =  sum_z2 + x_cm_value*2/10;
}

sum_y1 =  sum_y1 + sum_z1*8/10;
sum_y2 =  sum_y2 + sum_z2*8/10;
}

sum1 = sum1 + sum_y1*(4)/10;
sum2 = sum2 + sum_y2*(4)/10;
}

printf("first integral = %f, x_com = %f", sum1, sum2/sum1);
return 0;
}


In this, the entire domain is split into ranges of regular integrals according to the limits and by adding all the strips of z, one strip of y was obtained, similarly, by adding all the strips of y, one part of x is obtained which is summed up to obtain the final integral. Am I going wrong somewhere?

Probably classic example of Integer Division, maybe try adding dots . after integer numbers in the code to get the correct / intended fractional values:

#include <stdio.h>
#include <math.h>

double func(double x, double y, double z){
double out;
double condition;
//defining the domain
condition = z*z + pow((sqrt(x*x + y*y) - 3.), 2); // "." here not strictly necessary, but no harm done & indicating the double

if (condition <= 1. && x >= 0. && y >= -4.){// "." here not strictly necessary, but no harm done & indicating the double
out = 1.; // "." here not strictly necessary, but no harm done & indicating the double
}
else{
out = 0.; // "." here not strictly necessary, but no harm done & indicating the double
}

return out;
}

int main(){

int i,j,k;

double sum1 = 0, sum2 = 0;

for (i = 0; i < 10; i++){
double x = 0 + i*(4.)/10; // Note "." after 4
double sum_y1 =0, sum_y2 =0;

for (j=0;j<10;j++){
double y = -4 + j*(8.)/10; // Note "." after 8
double sum_z1 =0, sum_z2 =0;

for (k=0;k<10;k++){
double z = 1 - k*(2.)/10; // Note "." after 2
double x_cm_value = x*func(x,y,z);

sum_z1 =  sum_z1 + func(x,y,z)*2./10; // Note "." after 2
sum_z2 =  sum_z2 + x_cm_value*2./10; // Note "." after 2
}

sum_y1 =  sum_y1 + sum_z1*8./10; // Note "." after 8
sum_y2 =  sum_y2 + sum_z2*8./10; // Note "." after 8
}

sum1 = sum1 + sum_y1*(4.)/10; // Note "." after 4
sum2 = sum2 + sum_y2*(4.)/10; // Note "." after 4
}

printf("first integral = %f, x_com = %f", sum1, sum2/sum1);
return 0;
}