Physics C++

1. 0  This is a factorial code by which you can calculate the value of the factorial of 160 or more. This is really important  for the student of the Department of Physics......

#include <iostream>

#include<stdio.h>

#include<cmath>

using namespace std;

double fac( int n){              // The graet factorial finction

if (n==0) return 1;

else return n*fac(n-1);

}

int main ()

{

int i;

for(i=0; i<=100; i++){ 

cout << i <<" : "<<  fac(i)  << endl;

 

}

  return 0;

}

//////////////   For Double Factorial Function the code can be modified as..................


#include<iostream>
#include<cmath>


using namespace std;


double  fac( double n ){
      if (n==0) return 1;
      else return n*fac(n-1);
      }


double dfac( double n){
double val;
val=fac(n);
return fac(val);
      }
      
int main(){   
cout<<dfac(4)<<endl;
return 0;
}

 2.0  This is trick of Fourier, You may get this in your Mathematical Physics Book written by Arfken. If you have some graphing utility like gnuplot, you can get visible interpretation of this....

#include<iostream>
#include<stdio.h>
#include<cmath>
#include<stdlib.h>

using namespace std;
FILE*cplab;

double sinc(double x, int m){
double val=0.0;
int n=0;
for(n=1 ;n<=m;n++){
val+=pow(-1,n+1)*sin(n*x)/n;
}
return 2*val;
}

int main(){
cplab=fopen("gnu21","w");
double n=0;
double x=0.0, val=0.0,sgn=1.0;

for(x=0; x<=M_PI; x+=M_PI/80){

fprintf(  cplab, "%lf   %lf   %lf  %lf  %lf  %lf  %lf  %lf  \n",x,sinc(x,2),sinc(x,3),sinc(x,4),sinc(x,5),sinc(x,6),sinc(x,7),sinc(x,8));
}

fclose(cplab);
return 0;
}

3.0  This is trick for the visual interpretation of Legendre Function , You may get this in your Mathematical Physics Book written by Arfken. If you have some graphing utility like gnuplot, you can plot this.....

#include<iostream>
#include<cmath>

using namespace std;
FILE*cplab;

double legendre(double x, double n){                  // legendre function
if(n==0) return 1;
if(n==1) return x;
else return (2*x-x/n)*legendre(x,(n-1))-(1-1/n)*legendre(x,(n-2));
}

 

int main(){
cplab=fopen("lege2","w");

double i =0.0,  x=0.0;

for( x =-1; x<=1; x+=0.01){          // this x value is for graphing purpose

fprintf(cplab,"%lf \t %lf  \t  %lf  \t  %lf  \t  %lf  \t  %lf \n ", x , legendre(x,1), legendre(x,2),  legendre(x,3),  legendre(x,4), legendre(x,5) );
}

fclose(cplab);
return 0;
}

4.0  This is the code for calculating the square root of a positive real number.

#include<iostream>
#include<cmath>

using namespace std;

FILE*cplab;                               

double mysqroot( double b, int n){          //  square root calculator
  int i=0;
  double xn=1.0;                                     // xn is the approximated root        
  double val=0.0 , boo=0.0;   
  for(i=1; i<n; i++){
   
    val=0.5*(xn+b/xn);                            // Eq. for finding root
    xn=val; 
    boo=val-sqrt(b);                                // boo is the error between original and calculated value
   fprintf(cplab,"%d %lf \n",i,boo);         
    }
return val;
} 

int main(){
cplab=fopen("gnu21","w");               

int iter=30;                               // iteration is declared as 30
double b=0.0;

cout<<" Put the value:";
cin>>b;

printf("  My square root is:  %14.12lf \n",  mysqroot(b,iter));
printf("  Square root  from built in function :  %14.12lf \n", sqrt(b) );

fclose(cplab);                        

 return 0;
}

  

5.0  This function is like 

You will have find the value for n terms..

#include < iostream >
#include< cmath >

using namespace std;
       
       double mypi(int x , int n){
               if(n==1) return sqrt(x);
               else return  sqrt(2+ mypi( x , n-1));
               }            
int main ( )   { 
    int x=2,n=10;    // Here is the input value
   
    double val=0.0,nval=0.0, fval=0.0;
   
    val = powf(x,n);
    nval = mypi(x,n-1);
    fval = val*sqrt(x-nval);
   
    cout<< fval <<endl;
   
    return 0;
}





6.0 Calculate the square root of an equation using Bisection method. For theory you can look through any numerical book. Here is the direct code of Bisection method.....


#include<iostream>
#include<stdio.h>
#include<cmath>

using namespace std;

#define MAX 50                   // Maximum iteration number
#define EPS 1.0e-10            //define accuracy of convergence


double f(double x){ return (x+cos(x)); }   // f is the function whose argument is double

double bisect(double a,double b){
    int n;
    double c;
   
    for(n=0;n<MAX;n++){
        c=(a+b)/2.0;
        if(fabs(f(c))<EPS)    break;        //break if solution becomes close to the convergence
        if(f(a)*f(c)<0.0)    b=c;
        else a=c;
    }
    return c;
}

int main(void){
    double a,b,root;
   
    cout<<" Enter 2 starting values :  "<<endl;

    cout<<"A : ";
    cin>>a;
    cout<<"B : ";
    cin>>b;
    root = bisect(a,b);
   
    if (f(a)*f(b)>0.0) cout<<" the root is not located in this range "<<endl;
   
    else cout<<"the approximate Root is : "<<root<<endl;
   
    return 0;
}

7.0  Find the area bounded by two square of radius 1.0

Gnuplot drawing of the area

Here is the direct code,
#include<iostream>
#include<stdio.h>
#include<math.h>
#include<stdlib.h>

using namespace std;

FILE*cplab;                   

double getme(int N){
int i=0, seed=12345 , count=0;
double x=0.0 , y=0.0;

seed=time(NULL);
srand(seed);

for(i=1; i<=N ; i++){
x=(double)rand()/(double)RAND_MAX;
y=(double)rand()/(double)RAND_MAX;

if(    sqrt((x-1.0)*(x-1.0)+(y-0)*(y-0))<=1.0  
&&  sqrt((x-0)*(x-0)+(y-1.0)*(y-1.0))<=1.0
)
{
count++;
// fprintf(cplab,"%lf %lf \n",x,y);
}

}
return (double)count/(double)N;
}



int main(){

cplab=fopen("gnu3","w");             

int i=0 ;
double val=0.0;
long double ori=0.570796;

for(i=1000; i<=1000000; i*=10){

val=getme(i);

fprintf(cplab,"%d \t %10.8e \n",i,fabs(val-ori));  

printf("i=%d val=%14.12lf abs(val-0.570796)=%e \n",i,val,fabs(val-ori));
}

fclose(cplab);        
return 0;
}




8.0  Calculate the root of an equation using Newton Raphson method. For theory you can look through any numerical book. Here is the direct code of Newton Raphson method.....

#include<iostream>
#include<stdio.h>
#include<cmath>

using namespace std;

#define MAX 50                  // Maximum iteration number
#define EPS 1.0e-10            //define accuracy of convergence


double f(double x){ return (x*x*x-3*x-5); }  
 
double fdr(double x){ return ( 3*x*x-3 ); }

double newton(double a,double b){
    int n;
    double c;
    double val;
    c=(a+b)/2.0;
    
    for(n=0;n<MAX;n++){    
    val=c-(f(c)/fdr(c));
    c=val;    
    }
    return c;
}

int main(void){
    double a,b,root;
    
    cout<<"Enter 2 starting values : "<<endl;

    cout<<"A : ";
    cin>>a;
    cout<<"B : ";
    cin>>b;
    root = newton(a,b);
    
    if (f(a)*f(b)>0.0) cout<<" the root is not located in this range "<<endl;
    
    else printf("The approximate root is:%14.12lf\n",root);
  
    return 0;
}


9.0 Suppose You need to find the Determinant of a 4D square matrix.
Use Levi civita function that will give you the proper sign, and find the determinant.

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

using namespace std;



int levi(int i, int j, int k, int l)
{
    if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) return 0;
    else return ( (i-j)*(i-k)*(i-l)*(j-k)*(j-l)*(k-l)/12 );
}

int main(){
int val=0, i , j , k, l;
int AA[4][4]={ {1,2,3,4},{3,4,2,1},{2,5,6,8},{4,5,7,9} };

for(i=1 ; i<=4 ; i++){
for(j=1 ; j<=4; j++){
for(k=1 ; k<=4; k++){
for(l=1 ; l<=4; l++){
       
val=val+levi(i , j , k ,l)*AA[0][i-1]*AA[1][j-1]*AA[2][k-1]*AA[3][l-1];

//cout<<i<<j<<k<<l<<" : "<<levi(i,j,k,l)<<endl;
}
}
}
}
cout<<val<<endl;

return 0;

} //here we have divided that part by 12, because of 1/(1!*2!*3!) . for more information Look.

If you wish to find the determinant from online, click here

10.0 For a 5D square matrix. Find the determinant.

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

using namespace std;



int levi(int i, int j, int k, int l ,int m)
{
    if((i==j)||(i==k)||(i==l)||(i==m)||(j==k)||(j==l)||(j==m)||(k==l)||(k==m)||(l==m)) return 0;
    else return ( (i-j)*(i-k)*(i-l)*(i-m)*(j-k)*(j-l)*(j-m)*(k-l)*(k-m)*(l-m)/288 );
}


int main(){
int val=0, i , j , k, l, m;
int AA[5][5]={ {1,2,3,4,5},
               {3,4,2,1,6},
               {2,5,6,8,7},
               {4,5,7,9,8},
               {5,6,7,8,9} };

for(i=1 ; i<=5 ; i++){
for(j=1 ; j<=5; j++){
for(k=1 ; k<=5; k++){
for(l=1 ; l<=5; l++){
for(m=1 ; m<=5; m++){  

     
val=val+levi(i , j , k ,l , m)*AA[0][i-1]*AA[1][j-1]*AA[2][k-1]*AA[3][l-1]*AA[4][m-1];

//cout<<i<<j<<k<<l<<m <<" : "<<levi(i,j,k,l,m)<<endl;}
}
}
}
}

cout<<val<<endl;

return 0;

}

// here we have divided by 1!*2!*3!*4!. For more info look in Wiki.

11.0  Suppose a radioactive material is undergoing radioactive decay. The count per minute is between   600 ~800 .Show that the decay scheme follow Gaussian distribution.

Gnuplot drawing of  Random Nature

#include<iostream>
#include<stdio.h>
#include<cmath>
#include<stdlib.h>
#define pi 0.398942

using namespace std;
 FILE * cplab;

int main()                        // Expat No- 06, Random Nature of Radioactive decay
{
cplab = fopen("gnudata66", "w");
             
int i=0, j=0, k=0, seed=12345 , count=0,  N=280;     // Observation No
double x=0.0 , y=0.0, area=0.0 ;
double a=200.0, val=0.0, val2=0.0, xbar=0.0, sigma=0.0, Pb=0.0;

double arr[280], ab[280], absq[280],  tt[280];

seed = time(NULL);
srand( seed );

for( i=0;  i<N;  i++ )  {
x=a*(double)rand() / (double)RAND_MAX;
arr[i]=600.0+x;                                       // Count rate is between 600~800

val=val+arr[i];
}

xbar = val/280;

for( j=0;  j<N;  j++ ){
         ab[ j ] = arr[ j ] - xbar;
         absq[ j ] = ab[j]*ab[ j ];
         val2 = val2+absq[ j ];
         }
         sigma = sqrt(val2/279.0);

for(k=0; k<N; k++){

         tt[k] = ab[k] / sigma;

         Pb = pi*exp(-( tt[ k ]*tt[ k ] ) / 2.0);

         fprintf( cplab, "%lf  %lf  \n", tt[k], Pb);
         }
        
fclose(cplab);

return 0;
}

11.0 Suppose there are 9 integer value given. Arrange those integer in Ascending Order.

#include<iostream>
#include<stdio.h>
#include<stdlib.h>

using namespace std;

int main(){
  int arr[9]={-2,4,5,-7,8,9,12,-13,27};
  int x=0, tmp=0  , indx=0, i=0,j=0;
 

  for(j=0; j<9; j++)
    {
      indx=0;
      x=arr[j];
     
      for(i=j+1; i<9; i++)
    {
      if(x>arr[i]){    x=arr[i];    indx=i; }
    }
     
     
      //printf(" minimum value =%d \n",x);
      //printf(" minimum value  position =%d \n",indx); 
     

if(indx !=0)
    {
      tmp=arr[j];
      arr[j]=arr[indx];
      arr[indx]=tmp;
    }
     
    }
  cout<<"Numbers in Ascending order: ";
  for(i=0; i<9; i++) cout<<arr[i]<<" ";
  cout<<endl;

return 0;
}