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

double f(double z,double M,double Cc,double R){
 return pow(R*pow((1+z),4)+M*pow((1+z),3)+Cc+(1-R-M-Cc)*pow((1+z),2),-0.5);}

double trapezio(float z,double M,double Cc,double R){
 double Tk=0, T0, dx, soma, p=1.e-6;
 int j=0, i;
 do{
  dx=(z-0)/pow(2,j);
  if(j==0) T0=0.5*(z-0)*(f(0,M,Cc,R)+f(z,M,Cc,R));
  else T0=Tk;
  soma=0;
  for(i=1;i<=pow(2,j)-1;i+=2) soma+=f(0+i*dx,M,Cc,R);
  Tk=T0/2+dx*soma;
  j++;
  }
 while(fabs(Tk-T0)>p);
 return Tk;
 }

int main(){
 FILE *fentrada = NULL;
 fentrada = fopen&#40;"Distanciacomovel.dat","w"&#41;; 
 double M1,M2,M3,Cc1,Cc2,Cc3,R1,R2,R3,z;
 printf("Este programa calcula a idade do universo dados os parametros de densidade e o redshift desejado\n"); 
 printf("Entre com os valores dos parametros de densidade do modelo 1(M / Cc / R)\n");
 scanf("%lf\n %lf\n %lf\n",&M1;,&Cc1;,&R1;);
 printf("Entre com os valores dos parametros de densidade do modelo 2(M / Cc / R)\n");
 scanf("%lf\n %lf\n %lf\n",&M2;,&Cc2;,&R2;);
 printf("Entre com os valores dos parametros de densidade do modelo 3(M / Cc / R)\n");
 scanf("%lf\n %lf\n %lf\n",&M3;,&Cc3;,&R3;);
 double zlimite =50 , passo=0.5, n=0 , s;
 for(z=0;z<zlimite; z+=passo){
  fprintf(fentrada, "%lf\t %lf\t %lf\t %lf\n",z ,trapezio(z,M1,Cc1,R1),trapezio(z,M2,Cc2,R2),trapezio(z,M3,Cc3,R3));
  s=n*100/(zlimite*(1/passo));
  printf("%.2lf\n" , s);
  n+=1;
 }

 fclose(fentrada);
 return 0; 
}