<分区>
我正在编写一个 Runge-Kutta ODE 求解器，它接受初始条件和“事件”并求解方程。它在这些事件处开始和停止以形成“补丁”。然后将这些“补丁”放入保存的数组中。
该程序在对大量步骤进行评估时运行良好，但是，当对少量步骤进行评估时，它给出了奇怪的结果，该结果会在程序评估之间发生变化，有时会非常大(~e200) 或极小 (~e-200)。
我试图用 printf 语句来探究这个问题，似乎变量的值在调用 print 语句时发生了变化。
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] 在此引用之前为其分配了内存。
我在将数组解析为函数时使用了很多 malloc，我认为这可能是错误发生的地方。我最初试图尽量减少对 malloc 的调用，因为我知道它很慢，所以我正在解析指向函数的指针，但是，我现在已经换成在函数中尽可能多地声明(这就是为什么有些东西被声明了两次，但只有一个应该在范围内)用于调试。
我认为这可能与 6.0 的类型转换有关，并且它可能被零除，但事实并非如此，因为输出有时也非常小。它不稳定的事实似乎符合我的想法，即它是某种内存泄漏问题。我试着用 gdb 调试它，但我无法让它在我的电脑上运行。
我在运行 OS 10.13.6 的 Mac 上用 gcc 8.2 编译
如果有人有任何想法可能是错误的，我将非常感激。我不知道其他人是否遇到过类似的问题，我找不到任何其他具有相同问题的问题，但后来我不知道该寻找什么。非常感谢。
此处引用所有代码:
#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;
}