I'm writing a Runge-Kutta ODE solver which takes in initial conditions, and 'events' and solves the equation. It starts and stops at these events to form 'patches'. These 'patches' are then put into an array which is saved.
The programme worked fine when it was evaluated for a large number of steps, however, when it was evaluated for a small number of steps it gave strange results, which would change in value between evaluations of the programme, sometimes being extremely large (~e200) or extremely small (~e-200).
I tried to probe the issue with printf statements, and it seems that the value of the variables changes on calling the print statement.
printf("after solution K1_M[0] = %e,K1_M[1] = %e,K1_M[2] = %e\n",k1_M[0],k1_M[1],k1_M[2]);//This is fine
printf("Going into K2 Mx: %e, My: %e, Mz: %e\n",M_f[0],M_f[1],M_f[2]);// This is fine
printf("What is happening to this division? k1_M[0] = %e, 1/6.0 = %e, k1_M[0]/6.0 = %e\n",k1_M[0],1/6.0,k1_M[0]/6.0);
int hh = 1,hhh=3;
hhh = hh + hhh;
printf("Print statement to make sure other statements exicute in order %d\n",hhh);
M_f[0] = M_inf[1] + k1_M[0]/(double)(6.0);//This is not fine
printf("What is happening to this division? k1_M[0] = %e, 1/6.0 = %e, k1_M[0]/6.0 = %e\n",k1_M[0],1/6.0,k1_M[0]/6.0);
Example of the execution of the code and issue present
k1_M[0] has memory allocated to it just before this quotation.
I've used malloc a lot when parsing arrays into the functions, and I think that this is probably where the error occurs. I was initially trying to minimise the calls to malloc, as I understand that it is quite slow, so I was parsing the pointers to the function, however, I have now swapped to declaring as much as possible in the function (which is why some things are declared twice, but only one should be in scope) for debugging.
I thought that it could be something to do with the typecasting of 6.0, and that it might be dividing by zero, but that can't be it as the output is sometimes very small too. The fact that it is inconstant seems to go with my idea of it being some kind of memory leakage issue. I tried to debug it with gdb but I couldn't get it to run on my computer.
I was compiling with gcc 8.2, on a mac running OS 10.13.6
If anyone has any ideas what might be wrong I'd be extremely grateful. I don't know if other people have had similar problems, I couldn't find any other questions with the same issue, but then I didn't know exactly what to look for. Thanks very much.
All of the code is quoted here:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#include "plcdemos.h"
//#include <unistd.h>
//#include <string.h>
double vec_max(double *vec, int length){//find maximum value in a double vector of length length
double max_so_far = vec[0];
for(int i = 1;i<length;i++){
if(vec[i] > max_so_far){
max_so_far = vec[i];
}
}
return max_so_far;
}
double *bloch_patch(double t, double M[3], double Amplitude_RF, double gyro, double T_1, double T_2, double M_0, double B_z,int pulse_on){
//(M_inf[0],M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on)
double B_x = 0.0;
double w_l = B_z*gyro;
double w = w_l;
double w_1 = 0.0;
double dM[3] = {0.0,0.0,0.0};
double (*dM_pointer);// malloc(3*sizeof(double));
dM_pointer = (double*) malloc(3*sizeof(double));
if(pulse_on == 1)
B_x = Amplitude_RF;
//printf("Value of B_x in function is %e\n",B_x);
w_1 = gyro*B_x;
//w_l = B_z*gyro;
//ODE
//printf("Mx: %e, My: %e, Mz: %e\n",M[0],M[1],M[2]);
//printf("dM[1] First term is: %e\n",(w_l - w)*M[1]);
//printf("dM[1] Second term is: %e\n",- M[0]/T_2);
//printf("dM[2] First term is: %e\n",-(w_l - w)*M[0]);
//printf("dM[2] Second term is: %e\n",w_1*M[2]);
//printf("dM[2] Third term is: %e\n",- M[1]/T_2);
//actual ODE
dM[0] = (w_l - w)*M[1] - M[0]/T_2;
dM[1] = -(w_l - w)*M[0] + w_1*M[2] - M[1]/T_2;
dM[2] = -w_1*M[1] - (M[2] - M_0)/T_1;
//printf("dM values are: %e, %e, %e\n",dM[0],dM[1],dM[2]);
//end of ODE
//http://www.mmrrcc.upenn.edu/mediawiki/images/d/d1/Bloch.pdf
//printf("words %lf\n",dM[0]);
dM_pointer = dM;
return dM_pointer;
//return dM[0];
}
double *JVerner_step(double actual_dt,double *a,double *b,double *c,double *d,double *e,double *f,double *g, double M_inf[4],double Amplitude_RF,double gyro,double T_1,double T_2,double M_0,double B_z,int pulse_on){
//double *JVerner_step(double actual_dt,double *k1_M,double *k2_M,double *k3_M,double *k4_M,double *k5_M,double *k7_M,double *k8_M, double M_inf[4],double Amplitude_RF,double gyro,double T_1,double T_2,double M_0,double B_z,int pulse_on){
//printf("Jstep6\n");
double *M_out;
M_out = (double*)malloc(3*sizeof(double));
double M_f[3] = {0.0,0.0,0.0};
//was previously parsing kn_Ms into function. Where a-g are passed in. a-g do nothing.
double *k1_M,*k2_M,*k3_M,*k4_M,*k5_M,*k7_M,*k8_M;
k1_M = (double*)malloc(3 * sizeof(double)); //Try: does allocating the memory in the function make it work?
k2_M = (double*)malloc(3 * sizeof(double)); //It doens't fix it, but it doesn't make it less clear either
k3_M = (double*)malloc(3* sizeof(double));
k4_M = (double*)malloc(3 * sizeof(double));
k5_M = (double*)malloc(3* sizeof(double));
//k6_M not needed
k7_M = (double*)malloc(3 * sizeof(double));
k8_M = (double*)malloc(3 * sizeof(double));
for(int j = 0;j<3;j++){
k1_M[j] = 0.0;//force set to be zero
k2_M[j] = 0.0;// not using these
k3_M[j] = 0.0;
k4_M[j] = 0.0;
k5_M[j] = 0.0;
k7_M[j] = 0.0;
k8_M[j] = 0.0;
}
printf("Actual_dt: %e\n",actual_dt);
printf("Mx %e, My %e, Mz %e\n",M_inf[1],M_inf[2],M_inf[3]);
//M_inf[0] = M_t[2*i + 1];//time
//M_inf[1] = M_t[2*i + 2];//space
printf("before solution K1_y = %e\n",k1_M[1]);
M_f[0] = M_inf[1];
M_f[1] = M_inf[2];
M_f[2] = M_inf[3];
k1_M = bloch_patch(M_inf[0],M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k1_M[0] = k1_M[0]*actual_dt;//done
k1_M[1] = k1_M[1]*actual_dt;
k1_M[2] = k1_M[2]*actual_dt;
printf("after solution K1_M[0] = %e,K1_M[1] = %e,K1_M[2] = %e\n",k1_M[0],k1_M[1],k1_M[2]);//This is fine
printf("Going into K2 Mx: %e, My: %e, Mz: %e\n",M_f[0],M_f[1],M_f[2]);// This is fine
printf("What is happening to this division? k1_M[0] = %e, 1/6.0 = %e, k1_M[0]/6.0 = %e\n",k1_M[0],1/6.0,k1_M[0]/6.0);
int hh = 1,hhh=3;
hhh = hh + hhh;
printf("Print statement to make sure other statements exicute in order %d\n",hhh);
M_f[0] = M_inf[1] + k1_M[0]/(double)(6.0);//This is not fine
printf("What is happening to this division? k1_M[0] = %e, 1/6.0 = %e, k1_M[0]/6.0 = %e\n",k1_M[0],1/6.0,k1_M[0]/6.0);
M_f[1] = M_inf[2] + k1_M[1]/(double)(6.0);
M_f[2] = M_inf[3] + k1_M[2]/(double)(6.0);
printf("Going into K2 Mx: %e, My: %e, Mz: %e\n",M_f[0],M_f[1],M_f[2]);
k2_M = bloch_patch(M_inf[0]+(actual_dt/6.0),M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k2_M[0] = k2_M[0]*actual_dt;//done
k2_M[1] = k2_M[1]*actual_dt;
k2_M[2] = k2_M[2]*actual_dt;
printf("K2_y = %e\n",k2_M[1]);
M_f[0] = M_inf[1] + (4.0*k1_M[0] + 16.0*k2_M[0])/75.0;
M_f[1] = M_inf[2] + (4.0*k1_M[1] + 16.0*k2_M[1])/75.0;
M_f[2] = M_inf[3] + (4.0*k1_M[2] + 16.0*k2_M[2])/75.0;
k3_M = bloch_patch(M_inf[0]+(actual_dt*4.0/15.0),M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k3_M[0] = k3_M[0]*actual_dt;//done
k3_M[1] = k3_M[1]*actual_dt;// coeffs equal to 0.266666666666667
k3_M[2] = k3_M[2]*actual_dt;
printf("K3_y = %e\n",k3_M[1]);
M_f[0] = M_inf[1] + (5.0*k1_M[0]/6.0 - 8.0*k2_M[0]/3.0 + 5.0*k3_M[0]/2.0);
M_f[1] = M_inf[2] + (5.0*k1_M[1]/6.0 - 8.0*k2_M[1]/3.0 + 5.0*k3_M[1]/2.0);
M_f[2] = M_inf[3] + (5.0*k1_M[2]/6.0 - 8.0*k2_M[2]/3.0 + 5.0*k3_M[2]/2.0);
k4_M = bloch_patch(M_inf[0]+(actual_dt*2.0/3.0),M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k4_M[0] = k4_M[0]*actual_dt;//done
k4_M[1] = k4_M[1]*actual_dt;//coeffs equal to 0.666666666666667
k4_M[2] = k4_M[2]*actual_dt;
printf("K4_y = %e\n",k4_M[1]);
M_f[0] = M_inf[1] + (-165.0*k1_M[0]/64.0 + 55.0*k2_M[0]/6.0 - 425.0*k3_M[0]/64.0 + 85.0*k4_M[0]/96.0);
M_f[1] = M_inf[2] + (-165.0*k1_M[1]/64.0 + 55.0*k2_M[1]/6.0 - 425.0*k3_M[1]/64.0 + 85.0*k4_M[1]/96.0);
M_f[2] = M_inf[3] + (-165.0*k1_M[2]/64.0 + 55.0*k2_M[2]/6.0 - 425.0*k3_M[2]/64.0 + 85.0*k4_M[2]/96.0);
k5_M = bloch_patch(M_inf[0]+actual_dt*5.0/6.0,M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k5_M[0] = k5_M[0]*actual_dt;//done
k5_M[1] = k5_M[1]*actual_dt;//coeffs equal to 0.833333333333333
k5_M[2] = k5_M[2]*actual_dt;
printf("K5_y = %e\n",k5_M[1]);
//no need for k6
M_f[0] = M_inf[1] + (-8263.0*k1_M[0]/15000.0 + 124.0*k2_M[0]/75.0 - 643.0*k3_M[0]/680.0 - 81.0*k4_M[0]/250.0 + 2484.0*k5_M[0]/10625.0);
M_f[1] = M_inf[2] + (-8263.0*k1_M[1]/15000.0 + 124.0*k2_M[1]/75.0 - 643.0*k3_M[1]/680.0 - 81.0*k4_M[1]/250.0 + 2484.0*k5_M[1]/10625.0);
M_f[2] = M_inf[3] + (-8263.0*k1_M[2]/15000.0 + 124.0*k2_M[2]/75.0 - 643.0*k3_M[2]/680.0 - 81.0*k4_M[2]/250.0 + 2484.0*k5_M[2]/10625.0);
k7_M = bloch_patch(M_inf[0]+(actual_dt/15.0),M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k7_M[0] = k7_M[0]*actual_dt; //coeffs equal to 0.066666666666667
k7_M[1] = k7_M[1]*actual_dt;
k7_M[2] = k7_M[2]*actual_dt;
printf("K7_y = %e\n",k7_M[1]);
M_f[0] = M_inf[1] + (3501.0*k1_M[0]/1720.0 - 300.0*k2_M[0]/43.0 + 297275.0*k3_M[0]/52632.0 - 319.0*k4_M[0]/2322.0 + 24068.0*k5_M[0]/84065.0 + 3850.0*k7_M[0]/26703.0);
M_f[1] = M_inf[2] + (3501.0*k1_M[1]/1720.0 - 300.0*k2_M[1]/43.0 + 297275.0*k3_M[1]/52632.0 - 319.0*k4_M[1]/2322.0 + 24068.0*k5_M[1]/84065.0 + 3850.0*k7_M[1]/26703.0);
M_f[2] = M_inf[3] + (3501.0*k1_M[2]/1720.0 - 300.0*k2_M[2]/43.0 + 297275.0*k3_M[2]/52632.0 - 319.0*k4_M[2]/2322.0 + 24068.0*k5_M[2]/84065.0 + 3850.0*k7_M[2]/26703.0);
//M_f[0] = M_inf[1] + (3501.0*k1_M[0]/1720.0 - 300.0*k2_M[0]/43.0 + 297275.0*k3_M[0]/52632.0 - 319.0*k4_M[0]/2322.0 + 24068.0*k5_M[0]/84065.0 + k7_M[0]*0.144178556716475);
//M_f[1] = M_inf[2] + (3501.0*k1_M[1]/1720.0 - 300.0*k2_M[1]/43.0 + 297275.0*k3_M[1]/52632.0 - 319.0*k4_M[1]/2322.0 + 24068.0*k5_M[1]/84065.0 + k7_M[1]*0.144178556716475);
//M_f[2] = M_inf[3] + (3501.0*k1_M[2]/1720.0 - 300.0*k2_M[2]/43.0 + 297275.0*k3_M[2]/52632.0 - 319.0*k4_M[2]/2322.0 + 24068.0*k5_M[2]/84065.0 + k7_M[2]*0.144178556716475);
k8_M = bloch_patch(M_inf[0]+(actual_dt),M_f,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
k8_M[0] = k8_M[0]*actual_dt;// coeffs equal to 1.000 000 000 000 001 = 1 + 1e-15
k8_M[1] = k8_M[1]*actual_dt;
k8_M[2] = k8_M[2]*actual_dt;//0.144178556716475
printf("K8_y = %e\n",k8_M[1]);
// sixth order inc
M_out[0] = (3.0*k1_M[0]/40.0 + 875.0*k3_M[0]/2244.0 + 23.0*k4_M[0]/72.0 + 264.0*k5_M[0]/1955.0 + 125.0*k7_M[0]/11592.0 + 43.0*k8_M[0]/616.0);
M_out[1] = (3.0*k1_M[1]/40.0 + 875.0*k3_M[1]/2244.0 + 23.0*k4_M[1]/72.0 + 264.0*k5_M[1]/1955.0 + 125.0*k7_M[1]/11592.0 + 43.0*k8_M[1]/616.0);
M_out[2] = (3.0*k1_M[2]/40.0 + 875.0*k3_M[2]/2244.0 + 23.0*k4_M[2]/72.0 + 264.0*k5_M[2]/1955.0 + 125.0*k7_M[2]/11592.0 + 43.0*k8_M[2]/616.0);
printf("M_y = %e\n",M_out[1]);
// coeffs equal to 1.0 to machine precision
return M_out;
}
double *JsolveVernerPatch(double t,double t_start,double t_end, int n, double M_in[3],double Amplitude_RF, double gyro, double T_1, double T_2, double M_0, double B_z,int pulse_on){
//takes parameters of the function to be solved as well as the start point, the end point, and the regular stepsize
double actual_dt = (t_end-t_start)/(double)n;
printf("Actual_dt: %e\n",actual_dt);
double M_inf[4];//first for time, then three for mx my mz
double *k1_M,*k2_M,*k3_M,*k4_M,*k5_M,*k7_M,*k8_M; //k values in solver
double *F_out; //solution out
double *M_t;
M_t = (double*)malloc((4*(n+1)) * sizeof(double)); //vector will store points, t and M
//These are not currently used, how allocating inside of '_step' function
k1_M = (double*)malloc(3 * sizeof(double)); // allocation of memeory for k solver values I don't know
k2_M = (double*)malloc(3 * sizeof(double)); // if this is faster, but it seems like a sort of sensible thing
k3_M = (double*)malloc(3 * sizeof(double));
k4_M = (double*)malloc(3 * sizeof(double));
k5_M = (double*)malloc(3 * sizeof(double));
//k6_M not needed
k7_M = (double*)malloc(3 * sizeof(double));
k8_M = (double*)malloc(3 * sizeof(double));
F_out = (double*)malloc(3 * sizeof(double));
M_t[0] = t_start;
M_t[1] = M_in[0];
M_t[2] = M_in[1];
M_t[3] = M_in[2];
for(int i = 0;i<n;i++){
M_inf[0] = M_t[4*i + 0]; //time
M_inf[1] = M_t[4*i + 1]; //space
M_inf[2] = M_t[4*i + 2]; //space
M_inf[3] = M_t[4*i + 3]; //space
printf("Going into function: t: %e,Mx: %e,My: %e,Mz: %e\n",M_inf[0],M_inf[1],M_inf[2],M_inf[3]);
for(int j = 0;j<3;j++){
k1_M[j] = 0.0;//force set to be zero
k2_M[j] = 0.0;// not using these
k3_M[j] = 0.0;
k4_M[j] = 0.0;
k5_M[j] = 0.0;
k7_M[j] = 0.0;
k8_M[j] = 0.0;
}
F_out = JVerner_step(actual_dt,k1_M,k2_M,k3_M,k4_M,k5_M,k7_M,k8_M,M_inf,Amplitude_RF,gyro,T_1,T_2,M_0,B_z,pulse_on);
M_t[4*(i + 1) + 0] = M_inf[0] + actual_dt; //time
M_t[4*(i + 1) + 1] = M_inf[1] + F_out[0]; //space
M_t[4*(i + 1) + 2] = M_inf[2] + F_out[1]; //space
M_t[4*(i + 1) + 3] = M_inf[3] + F_out[2]; //space
}
return M_t;
}
double *JsolveVernerWhole(double t_start,double t_end, int n, double M_in[3],double Amplitude_RF, double gyro, double T_1, double T_2, double M_0, double B_z,int *pulse_v_pt,double *pulse_t_pt,int length_pulse){
//takes parameters of the function to be solved as well as the start point, the end point, and the regular stepsize
//double actual_dt = (t_end-t_star)/n;
//double *k1_M,*k2_M,*k3_M,*k4_M,*k5_M,*k7_M; //k values in solver
double *M_t,*M_patch;
double t;
int pulse_on = 0;
//typical_dt = dt;//pow(dt_min * dt_max,0.5);
int (*pulse_v);
double (*pulse_t);
pulse_v = (int*)malloc((length_pulse+1) * sizeof(int)); //how to get array of unknown length into function
pulse_t = (double*)malloc((length_pulse+1) * sizeof(double));
//pulse sequence
printf("Pulse length is %d\n",length_pulse);
for(int i=0;i<length_pulse+1;i++){
pulse_v[i] = pulse_v_pt[i]; //might be able to use use the pointers directly and skip this out
pulse_t[i] = pulse_t_pt[i];
printf("Pulse t is %e\n",pulse_t[i]);
//printf("pulse_v[%d] = %e\n",i,pulse_v[i]);
}
M_t = (double*)malloc((4*(n*(length_pulse)) + 4 ) * sizeof(double)); //vector will store points, t and M
printf("Function Whole M_t allocated size: %d\n",(4*(n*(length_pulse)) + 4 ));
M_patch = (double*)malloc(4*n* sizeof(double));
M_t[0] = t_start;
M_t[1] = M_in[0];
M_t[2] = M_in[1];
M_t[3] = M_in[2];
int start_at = 4;
printf("Time is: %e, Mx is: %e, My is: %e, Mz is: %e\n",M_t[0],M_t[1],M_t[2],M_t[3]);
for(int i = 0;i<(length_pulse);i++){
t = M_t[i*n*4 + 0];
//printf("Starting patch at t = %e\n",t);
M_in[0] = M_t[i*n*4 + 1];
M_in[1] = M_t[i*n*4 + 2];
M_in[2] = M_t[i*n*4 + 3];
if(length_pulse > 1){
for(int i = 0;i<length_pulse;i++){
if(t<=pulse_t[i]){
if(pulse_v[i] == 1){
pulse_on = 1;
}
break;
}
}
}//end of pulse sequence
printf("Tstart at %e, Tend at %e for current patch\n",pulse_t[i],pulse_t[i+1]);
printf("Index value is: %d\n",i);
M_patch = JsolveVernerPatch(t,pulse_t[i],pulse_t[i+1],n, M_in,Amplitude_RF, gyro, T_1, T_2, M_0, B_z, pulse_on);
printf("Time is: %e, Mx is: %e, My is: %e\n",M_patch[4*0 + 0],M_patch[4*0+ 1],M_patch[4*0 + 2]);
for(int j=1;j<n+1;j++){
printf("Time is: %e, Mx is: %e, My is: %e\n",M_patch[4*j + 0],M_patch[4*j + 1],M_patch[4*j + 2]);
/*
M_t[i*(n-1)*4 +4+ 0 + 4*j] = M_patch[4*j + 0];// time
M_t[i*(n-1)*4 +4+ 1 + 4*j] = M_patch[4*j + 1];
M_t[i*(n-1)*4 +4+ 2 + 4*j] = M_patch[4*j + 2];
M_t[i*(n-1)*4 +4+ 3 + 4*j] = M_patch[4*j + 3];// copy patch into matrix
*/
M_t[start_at + 0] = M_patch[4*j + 0];// time
M_t[start_at + 1] = M_patch[4*j + 1];
M_t[start_at + 2] = M_patch[4*j + 2];
M_t[start_at + 3] = M_patch[4*j + 3];// copy patch into matrix
printf("start_at: %d\n",start_at);
start_at = start_at + 4;
}
//for(int j=1;j<n+1;j++){
//printf("Time is: %e, Mz is: %e\n",M_patch[4*j + 0],M_patch[4*j + 3]);
/*
M_t[i*(n-1)*4 +4+ 0 + 4*j] = M_patch[4*j + 0];// time
M_t[i*(n-1)*4 +4+ 1 + 4*j] = M_patch[4*j + 1];
M_t[i*(n-1)*4 +4+ 2 + 4*j] = M_patch[4*j + 2];
M_t[i*(n-1)*4 +4+ 3 + 4*j] = M_patch[4*j + 3];// copy patch into matrix
*/
//}
}
/*
printf("Check loop\n");
for(int j=0;j<n+1;j++){
printf("Time is: %e, My is: %e\n",M_t[4*j + 0],M_t[4*j + 2]);
M_t[i*(n-1)*4 +4+ 0 + 4*j] = M_patch[4*j + 0];// time
M_t[i*(n-1)*4 +4+ 1 + 4*j] = M_patch[4*j + 1];
M_t[i*(n-1)*4 +4+ 2 + 4*j] = M_patch[4*j + 2];
M_t[i*(n-1)*4 +4+ 3 + 4*j] = M_patch[4*j + 3];// copy patch into matrix
}
*/
printf("start_at: %d\n",start_at);
return M_t;
}
int main(){
double Amplitude_RF,B_x,gyro,T_1,T_2,M_0,B_z,pulse_t[3];
int length_pulse;
int pulse_v[3];
double (*M_out);
Amplitude_RF = 1.0e-4;B_x = 0.0;gyro = 267.513e6;T_1 = 900e-3;T_2 = 90e-3;M_0 = 1.0;B_z = 3.5;
double M_in[3] = {0.0,1.0,0.0}; //xyz
double t_start = 0.0;
double t_end = 0.3;
int n = 3;
length_pulse = 1;
pulse_v[0] = 0; pulse_v[1] = 0; pulse_v[2] = 0;// pulse_v[1] = 0.0;
pulse_t[0] = t_start; pulse_t[1] = t_end; pulse_t[2] = t_end;
int *pulse_v_pt ;//= pulse_v;
double *pulse_t_pt;// = pulse_t;
pulse_v_pt = (int*)malloc((length_pulse + 1)*sizeof(int));
pulse_t_pt = (double*)malloc((length_pulse + 1)*sizeof(double));
//pulse sequence
printf("Pulse length is %d\n",length_pulse);
for(int i=0;i<length_pulse+1;i++){
pulse_v_pt[i] = pulse_v[i]; //might be able to use use the pointers directly and skip this out
pulse_t_pt[i] = pulse_t[i];
printf("Pulse t is %e\n",pulse_t[i]);
//printf("pulse_v[%d] = %e\n",i,pulse_v[i]);
}
//double tol = 4.6416e-7;
//double tol = 2.1544e-7;
//solver_out = Jsolve_RK6(t_start,t_end, dt_guess, M_in,Amplitude_RF, B_x, gyro, T_1, T_2, M_0, B_z,pulse_v,pulse_t, length_pulse);
//solver_out = Jsolve46(tol,rel_tol,t_start,t_end, dt_guess,dt_min,dt_max,M_in,Amplitude_RF, B_x, gyro, T_1, T_2, M_0, B_z,pulse_v,pulse_t,length_pulse);
//solver_out = Jsolve_RK4(t_start,t_end, dt_guess,M_in,a);
//M_out = JsolveFehlberg45(tol,1.0,t_start,t_end, dt_in,dt_min,dt_max,M_in,Amplitude_RF, B_x, gyro, T_1, T_2, M_0, B_z,pulse_v,pulse_t,length_pulse);
printf("Going into function t_start is: %e, coming out of function %e \n",t_start,t_end);
printf("Pulse_t[0] s: %e, pulse_t[1] %e \n",pulse_t[0],pulse_t[1]);
M_out = JsolveVernerWhole(t_start,t_end, n, M_in,Amplitude_RF, gyro, T_1, T_2, M_0, B_z,pulse_v_pt,pulse_t_pt,length_pulse);
//int skip = (int)fmax(ceil(1000/M_out[0]),1);
//double *mx,*my,*mz,*t,*modmy,*modmz,*diff,max_diff;
//mx = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//my = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//mz = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//t = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//modmy = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//modmz = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//diff = malloc(((length_pulse-1)*n + 1) * sizeof(double));
//max_diff = 0.0;
printf("end point = %d\n",((length_pulse)*n+1));
//printf("does this bit work\n");
//for(int i=0;i<((length_pulse)*n +1);i++){
//printf("index %d\n",i);
//printf("%lf %lf %lf %lf\n", solver_out[4*i+1],solver_out[4*i +2],solver_out[4*i +3],solver_out[4*i +4]);
//fprintf(fp,"%lf %lf %lf %lf \n", M_out[4*i+1],M_out[4*i +2],M_out[4*i+3],M_out[4*i +4]);
//printf("points = %d\n",points);
//printf("My[%d] = %e\n",points,M_out[4*i+2]);
// t[points] = M_out[4*i+0];
// printf("Time points %3.16e\n",t[points]);
// mx[points] = M_out[4*i+1];
// my[points] = M_out[4*i+2];
// printf("My[%d] = %e\n",points,my[points]);
// mz[points] = M_out[4*i+3];
//printf("t = %e,My data: %e\n",t[points],my[points]);
//printf("Index number is: %d\n",i);
//modmy[points] = exp(-t[points]/T_2);
//modmz[points] = 1-exp(-t[points]/T_1);
//diff[points] = fmax(fabs(my[points]-modmy[points]),fabs(mz[points]-modmz[points]));
//max_diff = fmax(max_diff,diff[points]);
// points++;
//printf("Time is: %e My is: %e\n",t[points],my[points]);
//fprintf(fp,"%3.16e %3.16e %3.16e %3.16e \n", M_out[4*i+1],M_out[4*i +2],M_out[4*i+3],M_out[4*i +4]);
//fprintf(fp,"the word data\n");
//}
printf("Solver_finished\n");
printf("Plotting\n");
printf("total points: %d\n",(int)(4*((length_pulse)*n+1)));
//FILE *fp = fopen(strcat(getcwd(cwd,sizeof(cwd)),"data.csv"),"W");
printf("opening file to write to\n");
FILE *fp = fopen("Bloch_RKR.csv","w");
if (fp == NULL)
{
printf("Error opening file!\n");
exit(1);
}
//printf("fp %f\n",fp);
//fp ;
printf("writing data to file\n");
//for(int i=0;i<(points);i++){
//printf("index %d\n",i);
//printf("%lf %lf %lf %lf\n", solver_out[4*i+1],solver_out[4*i +2],solver_out[4*i +3],solver_out[4*i +4]);
//fprintf(fp,"%lf %lf %lf %lf \n", M_out[4*i+1],M_out[4*i +2],M_out[4*i+3],M_out[4*i +4]);
// printf("%3.16e %3.16e %3.16e %3.16e\n", t[i],mx[i],my[i],mz[i]);
// fprintf(fp,"%3.16e %3.16e %3.16e %3.16e\n", t[i],mx[i],my[i],mz[i]);
//fprintf(fp,"the word data\n");
//}
for(int i=0;i<((length_pulse)*n +1);i++){
//printf("index %d\n",i);
//printf("%lf %lf %lf %lf\n", solver_out[4*i+1],solver_out[4*i +2],solver_out[4*i +3],solver_out[4*i +4]);
//fprintf(fp,"%lf %lf %lf %lf \n", M_out[4*i+1],M_out[4*i +2],M_out[4*i+3],M_out[4*i +4]);
//printf("points = %d\n",points);
//printf("My[%d] = %e\n",points,M_out[4*i+2]);
printf("%3.16e %3.16e %3.16e %3.16e\n", M_out[4*i+0],M_out[4*i+1],M_out[4*i+2],M_out[4*i+3]);
fprintf(fp,"%3.16e %3.16e %3.16e %3.16e\n", M_out[4*i+0],M_out[4*i+1],M_out[4*i+2],M_out[4*i+3]);
//t[points] = M_out[4*i+0];
//printf("Time points %3.16e\n",t[points]);
//mx[points] = M_out[4*i+1];
//my[points] = M_out[4*i+1];
// printf("My[%d] = %e\n",points,my[points]);
//mz[points] = M_out[4*i+3];
//printf("t = %e,My data: %e\n",t[points],my[points]);
//printf("Index number is: %d\n",i);
//modmy[points] = exp(-t[points]/T_2);
//modmz[points] = 1-exp(-t[points]/T_1);
//diff[points] = fmax(fabs(my[points]-modmy[points]),fabs(mz[points]-modmz[points]));
//max_diff = fmax(max_diff,diff[points]);
//points++;
//printf("Time is: %e My is: %e\n",t[points],my[points]);
//fprintf(fp,"%3.16e %3.16e %3.16e %3.16e \n", M_out[4*i+1],M_out[4*i +2],M_out[4*i+3],M_out[4*i +4]);
//fprintf(fp,"the word data\n");
}
fclose(fp);
printf("data written\n");
return 0;
}
