Confronto strutture

Ciao ragazzi!
devo fare un programma in C che confronti 2 funzioni (la seconda è contrassegnata da S). poichè queste funzioni appartengono a un codice più grande per il quale è necessario che si utilizzino delle strutture, ho pensato di mettere gli elementi della struttura dentro un Arrey per facilitarne la differenza.
quando lo eseguo..NEL TXT COMPAIONO BESTEMMIONI BIBLICI!
se qualcuno puo darci un occhio...posso ricambiare con fisica e matematica (sono un astrofisica!)!

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


// 1) definisco le costanti

const double mucrit=0.;
const double Asun=7.95495e+13;;
const double rhelio=0.004633333;
const double Kp0=0.05;
const double A=1.;	
const double Z=1.;
double T0=938.;	
double Pc=1015.;
double k0=0.000100;
const double c=3e8;	
double V0=2.6e-06;		
double Vmax=5.0e-6;
int	  VSWModel=0;
double VswAngl=0.8660;                      
const double Rhelio=100.;
const double omega=3.03008e-6;



// 2) definisco le strutture

typedef struct { //contiene le 3 derivate parziali di una cordinata
        double dr;             //derivata parziale risp r
        double dmu;            //derivata parziale risp mu
        double dphi;           //derivata parziale risp phi
} derivata; 

typedef struct { //contiene le 3 coordinate e le 3 derivate
        double r;
        double mu; //cos(teta)
        double phi;
        derivata Dr;           //derivata di r (oggetto con 3 dati)
        derivata Dmu;          //derivata di mu (oggetto con 3 dati)
        derivata Dphi;         //derivata di phi (oggetto con 3 dati)
} vect; 

typedef struct {  
        vect r;
        vect mu;
        vect phi;
} tensor; 






// 3) definisco le funzioni


double solarwind(double mu);
vect campomagnetico(double r,double mu);
vect campomagneticoS(double r,double mu);

int main()
{
FILE *fp = fopen("checkB.txt", "w");


 double T=1.;	/* energy of particle	MeV	*/
 double contmu;
 double contr;
for (contr=0.;contr<100.;contr++){
    double r=contr;
//fprintf(fp,"il vettore differenza dei campi magnetici a r = %f è\n",r);
               
               for (contmu=0.;contmu<=10.;contmu++){
          double mu=contmu/10.;
          fprintf(fp,"il vettore differenza dei campi magnetici a r = %f r mu = %f è\n",r,mu);

vect Bmag;
Bmag=campomagnetico(r,mu);
double B[3];
B[0]=Bmag.r;
B[1]=Bmag.mu;
B[2]=Bmag.phi;

vect BmagS;
BmagS=campomagneticoS(r,mu);     
double BS[3];
BS[0]=BmagS.r;
BS[1]=BmagS.mu;
BS[2]=BmagS.phi;

double deltaB[3];
int i;
    for(i=0; i<3; i++) {
        deltaB[i]= B[i]-BS[i];
        fprintf(fp,"mari %e; ste %e differenza %e \t",B[i],BS[i],deltaB[i]);
        fprintf(fp,"\n");   
        }      
}   
}
fclose(fp);

 return 0;
}


/////////////////////// vento solare ///////////////////
/*Solar wind description*/      
double solarwind(double mu){
 switch (VSWModel){
   case 0: return V0;
   case 1: if (fabs(mu)>VswAngl){return Vmax;} 
           else return (V0*(1+fabs(mu)));
   default:return V0;   
   }
}

// ////////////////////////// CAMPO MAGNETICO \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\

vect campomagnetico(double r,double mu){

vect Bmag;
double Vsw;
Vsw=solarwind(mu);

double PolSign;
double DelDirac;
  if (mu>=mucrit) PolSign=+1.; 
             else PolSign=-1.;
  if (mu==mucrit) DelDirac=1.;
              else DelDirac=0.;

Bmag.r = (Asun/r/r)*PolSign;
Bmag.Dr.dr = -2* (Asun/r/r/r)*PolSign;
Bmag.Dr.dmu = (Asun/r/r)*DelDirac;
Bmag.Dr.dphi = (Asun/r/r)* PolSign;
Bmag.mu = 0.;
Bmag.Dmu.dr = 0.;
Bmag.Dmu.dmu = 0.;
Bmag.Dmu.dphi = 0.;
Bmag.phi =(Asun/r/r)*PolSign*(-omega*(r-rhelio)*pow((1-mu*mu),1/2))/Vsw;
Bmag.Dphi.dr = (Asun/r/r/r/r)*PolSign*omega*pow((1-mu*mu),1/2)/Vsw*(-r*r+2*(r-rhelio)*r);
Bmag.Dphi.dmu = ((Asun/r/r)*PolSign*(omega*(r-rhelio)*pow((1-mu*mu),-1/2)*mu)/Vsw)+((Asun/r/r)*DelDirac*(-omega*(r-rhelio)*pow((1-mu*mu),1/2))/Vsw); 
Bmag.Dphi.dphi = 0.;

return Bmag;

}

//////////////// CAMPO MAGNETICO STE ////////////////////////////////

vect campomagneticoS(double r,double mu){
    vect Bmag;
  double V_SW;
  double PolSign; //emulate the heavside function for changing polarity on and behind the Neutral sheet
  double DelDirac;
  double Asun;
  double delta_m=0.;
  if (mu>=0) PolSign=+1.; // flat neutral sheet
             else PolSign=-1.;
  if (mu==0.) DelDirac=1.;
              else DelDirac=0.;
  V_SW=solarwind(mu);
  
  
  //magnetic Field vector

  Bmag.r=Asun/(r*r)*PolSign;
  Bmag.mu= Asun*delta_m/(r*rhelio*sqrt(1-mu*mu)) ;
  Bmag.phi=-Asun/(r*r)*(omega*(r-rhelio)*sqrt(1-mu*mu))/V_SW*PolSign;
  //derivative
  Bmag.Dr.dr=- 2.*Asun*PolSign/(r*r*r);
  Bmag.Dr.dmu=-2.*Asun*DelDirac/(r*r);
  Bmag.Dr.dphi=0.;
  
  Bmag.Dmu.dr=- Asun*delta_m/(r*r*rhelio*sqrt(1-mu*mu)) ; 
  Bmag.Dmu.dmu= Asun*delta_m*mu/(r*rhelio*pow(1-mu*mu,1.5));
  Bmag.Dmu.dphi=0.;
  
  Bmag.Dphi.dr=-Asun*(r-2*rhelio)*sqrt(1-mu*mu)*omega*(-PolSign)/(r*r*r)/V_SW;
  Bmag.Dphi.dmu=- Asun*(r-rhelio)*omega*( ((mu*mu-1.)*PolSign*V_SW) +V_SW*(-mu*PolSign+2.*(mu*mu-1)*DelDirac))/(r*r*sqrt(1-mu*mu)*V_SW*V_SW);
  Bmag.Dphi.dphi=0.;
  
  return Bmag;  

}