查看: 3575|回复: 1|关注: 0

[已答复] (挑战)关于bloch方程现有C++程序,谁能改成matlab的?

[复制链接]

新手

5 麦片

财富积分


050


0

主题

3

帖子

0

最佳答案
发表于 2011-4-16 16:04:51 | 显示全部楼层 |阅读模式
#include <stdio.h>
#include <math.h>
#define pi 3.1415929535898
#define GAMMA 4257*2*pi  /*  radian/gauss  */
grad_filestring[80];
rf_filestring[80];
float freq, rfx, rfy;
/***********************************************************/
/*          */
/* subroutine to define the Bloch Equation            */
/*               */
/*         Note: the array starts with 1     */
/*          */
/***********************************************************/
void derivs(x, y, dydx)
float y[], dydx[], x;
{
/** Note: the array starts with 1  **/
dydx[1] = (y[2]*freq - GAMMA*rfy*y[3]);      /* dMx/dt  */
dydx[2] = (GAMMA*rfx*y[3] - y[1]*freq);      /* dMy/dt  */
dydx[3] = GAMMA*(rfy*y[1] - rfx*y[2]);      /* dMz/dt  */
}
/***********************************************************/
/*          */
/* subroutine to perform 4th order Runge_Kutta        */
/*               */
/*       modified from rk4.c routine in numerical recipes. */
/*       ANSI and C++ stuff are removed.     */
/*          */
/***********************************************************/
#define NRANSI
#include "nrutil.h"
void rk4(y, dydx, n, x, h, yout, derivs)
float y[], dydx[], x, h, yout[];
void (*derivs)();
int n;
{
}
#undef NRANSI

main(argc, argv)
int argc;
char **argv;
{
float tsp;  /* sampling time in us */
float mx0, my0, mz0;    /* initial magnetizations in arbitary units */
float omegamax;         /* maximum frequency offset in Herz  */
float zmax;  /* maximum z in cm   */
int np_z;  /* number of points for z-array */
int np_omega;  /* number of points for omega array */
int np_g;  /* number of points for gradient/rf waveforms */
float delta_z;  /* slice offset in cm */
float g_scale;  /* gradient scaling factor */
int raster;  /* raster size of the output files */
float phi;  /* phase offset in radians for slice offset */
int nn;   /* index for magnetization output */
void rk4(), derivs(), free_vector();   /* subroutines */
float h, x, *y, *dydx, *yout;          /* variables for subroutines */
float *G, *B1, *B1x, *B1y, *z, *Mx, *My, *Mz, *omega, *kz;  /* arrays  */
register int j,i,m, k;
int n=3;      /* number of coupled ordinary
        differential equations.
        For Bloch, n = 3   */
char *malloc();
                        
FILE *Mx_file;     /* output file for Mx  */
FILE *My_file;     /* output file for My  */
FILE *Mz_file;     /* output file for Mz  */
FILE *z_file;      /* output file for z-coordinates  */
FILE *omega_file;  /* output file for omega-coordinates  */
FILE *grad_file;   /* gradient input file (gauss/cm)  */
FILE *rf_file;     /* rf input file (gauss)  */
FILE *kz_file;     /* kz output file */
FILE *fopen();

/* loading input parameters   */
if (argc != 14)
    {
            fprintf(stderr, "usage: bloch_spsp grad_file rf_file tsp(us) mx0 my0 mz0 omega_max(hz) np_omega z(cm) np_z np_g delta_z(cm) g_scale\n");
            exit(1);
           }

        if ((grad_file = fopen(argv[1], "r")) == NULL)
    {
                fprintf(stderr, "can't open %s\n", argv[1]);
                exit(1);
    }
if ((rf_file = fopen(argv[2], "r")) == NULL)
           {
                fprintf(stderr, "can't open %s\n", argv[2]);
                exit(1);
           }
         sscanf(argv[3],"%f", &tsp);
  sscanf(argv[4],"%f", &mx0);
  sscanf(argv[5],"%f", &my0);
  sscanf(argv[6],"%f", &mz0);
  sscanf(argv[7],"%f", &omegamax);
  sscanf(argv[8],"%d", &np_omega);
  sscanf(argv[9],"%f", &zmax);
  sscanf(argv[10],"%d", &np_z);
  sscanf(argv[11],"%d", &np_g);
  sscanf(argv[12],"%f", &delta_z);
  sscanf(argv[13],"%f", &g_scale);
  if((np_omega <2) || (np_z <2))
{
    printf("np_z and np_omega must be larger than 1 - program dies.\n");
    exit(1);
}
if (((Mx_file = fopen("mx_file", "w")) == NULL)  ||
     ((My_file = fopen("my_file", "w")) == NULL)  ||
     ((Mz_file = fopen("mz_file", "w")) == NULL) )
{
  fprintf(stderr, "can't open magnetization output files.\n");
                exit(1);
        }
if (((z_file = fopen("z_file", "w")) == NULL)  ||
     ((omega_file = fopen("omega_file", "w")) == NULL) ||
     ((kz_file = fopen("kz_file", "w")) == NULL) )
{
  fprintf(stderr, "can't open output z and omega files.\n");
                exit(1);
        }
/* set up arrays */
    y = vector(1, n);
    dydx = vector(1, n);
    yout = vector(1,n);

/* Memory allocation */
G = (float *)malloc(np_g*sizeof(float));
B1 = (float *)malloc(np_g*sizeof(float));
B1x = (float *)malloc(np_g*sizeof(float));
B1y = (float *)malloc(np_g*sizeof(float));
z = (float *)malloc(np_z*sizeof(float));
omega = (float *)malloc(np_omega*sizeof(float));
raster = np_omega*np_z;
Mx = (float *)malloc(raster*sizeof(float));
My = (float *)malloc(raster*sizeof(float));
Mz = (float *)malloc(raster*sizeof(float));
printf("memory allocated.\n");
/* Read in external wave forms ****/
printf("reading the input grad file ...\n");
    if(fread(G, sizeof(float), np_g, grad_file)!= np_g)
  {
  fprintf(stderr,"bad read of gradient file.\n");
  exit(1);
  }
printf("reading the input rf file ...\n");
    if(fread(B1, sizeof(float), np_g, rf_file)!= np_g)
  {
  fprintf(stderr,"bad read of RF file.\n");
  exit(1);
  }
/* parameter scaling */
tsp = tsp*0.000001;
        omegamax = omegamax*2*pi;

/*  initialize the spatial array */
    for (i=0; i<np_z; i++)
    z = -zmax+delta_z+(float)i*(zmax*2)/(float)(np_z-1);
/*  calculating the max phase - Note: the k-space is reversely defined. */
    phi = 0.0;
    for(i=0; i<np_g; i++)
      {
G = G*g_scale;
phi = phi + GAMMA*G*delta_z*tsp;
       }
/* generating the B1x and B1y arrays with slice offset */
    for(i=0; i<np_g; i++)
      {
phi = phi-GAMMA*G*delta_z*tsp;
B1x = B1*cos(phi);
B1y = B1*sin(phi);  /* Note: the way to generate the complex
      /* RF is proved in GE Note Book #2 */
      }
/* Bloch Simulation   */
  for (k=0; k<np_omega; k++)
   {
     omega[k] = -omegamax + (float)k*(omegamax*2)/(float)(np_omega-1);
     if (k%16 == 0)
     printf(" ... working on freq # %d out of %d ...\n", k, np_omega);
for (i=0; i<np_z; i++)
{
     nn = k*np_z+i;
     y[1] = mx0;
     y[2] = my0;
     y[3] = mz0;   /* initiation */
      
      freq = (omega[k] + GAMMA*G[0]*z);
     rfx = B1x[0];
     rfy = B1y[0];
     x= 0.0;
     derivs(x, y, dydx);  /* initiation */
     for(j=1; j<np_g; j++)
  {
  freq = (omega[k] + GAMMA*G[j]*z);
  rfx = B1x[j];
  rfy = B1y[j];     /* to be used for the next call of derivs */
  x = (float)(j-1)*tsp;
  rk4(y, dydx, n, x, tsp, yout, derivs);
  y[1] = yout[1];
  y[2] = yout[2];
  y[3] = yout[3];   /* reset initial values */
  derivs(x, y, dydx); /* reset initial values */  
  }

  Mx[nn] = yout[1];
  My[nn] = yout[2];
  Mz[nn] = yout[3]; /* grab the output */

}
   }
    free_vector(y, 1, n);
    free_vector(dydx, 1, n);
    free_vector(yout, 1, n);

/** writing the results **/
fprintf(stderr, "writing out the results ...\n");
        if ((fwrite(Mx,sizeof(float), raster, Mx_file)!= raster)  ||
            (fwrite(My,sizeof(float), raster, My_file)!= raster)  ||
            (fwrite(Mz,sizeof(float), raster, Mz_file)!= raster)  ||
     (fwrite(z,sizeof(float),  np_z, z_file)!= np_z)       ||
     (fwrite(omega,sizeof(float), np_omega, omega_file)!= np_omega) )
                {
                        fprintf(stderr, "write outfile returned < 0.\n");
                        exit(1);
                }

fclose(grad_file);
fclose(rf_file);
fclose(Mx_file);
fclose(My_file);
fclose(Mz_file);
fclose(z_file);
fclose(omega_file);
}
/***********************************************************/
/*          */
/* other subroutines              */
/*               */
/*         modified from nrutil.c in Numerical Recipes     */
/*          */
/***********************************************************/
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
float *vector(nl, nh)
long nl, nh;
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) printf("allocation failure in vector().\n");
return v-nl+NR_END;
}
void free_vector(v, nl, nh)
float *v;
{
free((FREE_ARG) (v+nl-NR_END));
}

Balac-Chupin.pdf

330.78 KB, 下载次数: 405

回复主题 已获打赏: 0 积分

举报

新手

21 麦片

财富积分


050


6

主题

29

帖子

0

最佳答案
发表于 6 天前 | 显示全部楼层
请问这段c++有注释吗,具体要怎么求解Bloch方程呢?
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表