/* SAR Calculator ************************************************************ */
/* Mike Weisdorf, 7/99 							       */
/* Calculates deposited RF energy and average RF power in VNMR pulse sequences */
/* See also, sarcalc.txt or local CMRR web site for more information.	       */
/* Version 1.2   */
/* Version 1.3  Lance DelaBarre 9_30_99  Eliminated error warnings by casting  */
/*              readgate, readfile, calcint, calmaxpwr as int instead of void. */
/*              Cast getenv as a char pointer to eliminate lint warning in getpwr */
/* Version 1.4	Lance DelaBarre 3_16_00  Had problems with repetative remote file system reads */
/*		changed the order of readin to try 3 times, not just three reads in a row like before. */

#include <stdio.h>
#include <math.h>
#define	sardecouple(target)									\
	if ( ((target == 's') || (target == 'S')) || ((target == 'a') || (target == 'A')) || ((target == 'y') || (target == 'Y')) ) { sarcalc("hard",dpwr,at,DD,1); }
#define POUND_LTR ('#')
#define saravgtime(tr) fprintf(stdout,"\nTotal energy in joules = %f\nAvg. power in Watts = %f\n",sarenergy,sarenergy/(tr))
#define sarcalc(shape,power,width,channel,ntimes)		\
	sprintf(parrname,#power); bpv=-17; ardx=0;		\
	if ((nar=getarray(parrname,parray))>1) {		\
	while (ardx<=nar) { cpv=parray[ardx]; 			\
		 if (bpv<cpv) bpv=cpv; ardx++; }		\
	fprintf(stdout,"\n%s is arrayed. SARcalc uses largest value: %f\n",parrname,bpv); } \
	else bpv=power;				\
	sprintf(warrname,#width); bwv=0; ardx=0;		\
	if ((nar=getarray(warrname,warray))>1) {		\
	while (ardx<=nar) { cwv=warray[ardx]; 			\
		 if (bwv<cwv) bwv=cwv; ardx++; }		\
	fprintf(stdout,"%s is arrayed. SARcalc uses largest value: %f\n",warrname,bwv); } \
	else bwv=width;				\
	sarenergy+=ntimes* getpwr(shape,bpv,bwv,channel)
#define sargetpwr(shape,power,width,channel,sumname)	\
	sprintf(parrname,#power); bpv=-17; ardx=0;		\
	if ((nar=getarray(parrname,parray))>1) {		\
	while (ardx<=nar) { cpv=parray[ardx]; 			\
		 if (bpv<cpv) bpv=cpv; ardx++; }		\
	fprintf(stdout,"\n%s is arrayed. SARcalc uses largest value: %f\n",parrname,bpv); } \
	else bpv=power;				\
	sprintf(warrname,#width); bwv=0; ardx=0;		\
	if ((nar=getarray(warrname,warray))>1) {		\
	while (ardx<=nar) { cwv=warray[ardx]; 			\
		 if (bwv<cwv) bwv=cwv; ardx++; }		\
	fprintf(stdout,"%s is arrayed. SARcalc uses largest value: %f\n",warrname,bwv); } \
	else bwv=width;				\
	sumname=getpwr(shape,bpv,bwv,channel)
int readgate();
int readfile();
int calcint();
int calmaxpwr();
const int DD=0, TT=1;
double sarenergy=0;
double parray[32], warray[32];
char sar[MAXSTR], parrname[16], warrname[16];
FILE *fp;
int n_ppts, xn, nar, ardx=0;
float retpwr, maxpwr, duration, cwv, cpv, bwv, bpv;
struct pp_t {
	float phase;
	float amp;
	float time;
	float gate;
	float ampsqr;
} seq[65536];
float getpwr(char ipath[400], float imaxpwr, float iduration, int ixn)
{
	char path[400];
	char tmplin[200];
	int gatechk, junkoutput, existtry;
	char *getenv();     /* This line eliminates error messages about uncast pointers from getenv Lance D 10/99*/
	maxpwr = imaxpwr;
	duration = iduration;	
	xn = ixn;
	junkoutput = calmaxpwr();
/*	fprintf(stderr, "maxpwr %f   duration %f\n", maxpwr, duration);	*/

	existtry = 1;
	while ((existtry <= 3) && (existtry != 0)) {
          sprintf(path, "%s/shapelib/%s.RF", getenv("vnmruser"), ipath);
	  fp=fopen(path, "r");
	  if (fp == NULL) {      /* Not in user shapelib dir*/
	    fclose(fp);
	    sprintf(path, "/vnmr/shapelib/%s.RF", ipath);
	    fp=fopen(path, "r");
	    if (fp == NULL) {    /* Not in /vnmr/shapelib directory */
		fclose(fp);
		existtry++;
	    } else existtry = 0;
	  } else existtry = 0;
	}
	if (existtry != 0) {
	  fprintf(stderr, "Cannot find file %s\n", ipath);
	  exit(1);
	} 
	fprintf(stderr, "%s\n", path); 				
	n_ppts = 0;
	while(!feof(fp))
	{
		if (NULL == fgets(tmplin, 200, fp))
		 break;
		if (POUND_LTR == *tmplin)
		/* printf("%s", tmplin) */;
		else
		 n_ppts++;
	}
	rewind(fp);
	while (POUND_LTR == fgetc(fp)) fgets(tmplin, 200, fp);
	fscanf(fp, "%*f%*f%*f");
        do {
	gatechk = fgetc(fp);
	switch (gatechk) {
                case 48:	/* 0: use gate */
		junkoutput = readgate();
		break;
		case 49:	/* 1: use gate */
                junkoutput = readgate();
		break;
                case 10:	/* \n: no gate */
		junkoutput = readfile();
		break;
		case 32:	/* space: check next */
		break;
                default:
                fprintf(stderr, "Bad format in file %s\n", path);
                exit(0);
        	}
	} while (gatechk==32);
        junkoutput = calcint();
	fprintf(stderr, "energy %f\n", retpwr);
	return retpwr;
}
int readgate()
{		
	extern struct pp_t seq[65536];
	int n = 0;
	char tmplin[200];
	rewind(fp);
	while(n < n_ppts) {
		if (POUND_LTR == fgetc(fp)) fgets(tmplin, 200, fp);
		else {
		fscanf(fp, " %f %f %f %f \n", &seq[n].phase, &seq[n].amp, &seq[n].time, &seq[n].gate);
                seq[n].ampsqr = seq[n].amp * seq[n].amp * seq[n].gate;
		n++; }
	}
	fclose(fp);
	return(1);
}
int readfile()
{	
	extern struct pp_t seq[65536];
	int n = 0;
	char tmplin[200];
 	rewind(fp);
 	while(n < n_ppts) {
		if (POUND_LTR == fgetc(fp)) fgets(tmplin, 200, fp);
		else {
		fscanf(fp, " %f %f %f \n", &seq[n].phase, &seq[n].amp, &seq[n].time);
		seq[n].gate = (float)1.00;
		seq[n].ampsqr = seq[n].amp * seq[n].amp;
		n++; }
	}
	fclose(fp);
	return(1);
}
int calcint()
{
	extern struct pp_t seq[65536];
	float ans = 0.0, rte = 0.0;
	int m = 0;
	float normw=0;
	float chkw;
	while(m < n_ppts) {
     		ans = ans+(seq[m].ampsqr * seq[m].time);
		rte = rte+seq[m].time;
		chkw = seq[m].ampsqr;
		if (normw < chkw) normw = chkw;
		m++;
	}
/*	fprintf(stderr, "calcint %f\n", ans); 				
	fprintf(stderr, "dur %f maxpwr %f\n", duration, maxpwr);
        fprintf(stderr, " normt %f normw %f\n", rte, normw);		*/
	retpwr = (ans * maxpwr * duration) / (normw * rte);
	return(1);
}
int calmaxpwr()
{
	float tcal[2];
	FILE *calfile;
	char calpath[100+MAXSTR], usn[MAXSTR];
	if (!xn) { getstr("dn", usn);
	 sprintf(calpath, "/vnmr/imaging/sartables/sarcalib.dn.%s", usn); }
	else { getstr("tn", usn); 
	 sprintf(calpath, "/vnmr/imaging/sartables/sarcalib.tn.%s", usn); }
	if ((calfile=fopen(calpath, "r"))==NULL) {
		fprintf(stderr, "Cannot find correct sarcalib for %s\n", usn);
		exit(1); }
/*	fprintf(stderr, "%s\n", calpath);				*/
	tcal[0]=-17;
	while ((tcal[0] < maxpwr) && !feof(calfile)) fscanf(calfile, " %f %f \n", &tcal[0], &tcal[1]);
	maxpwr = exp(0.2302585093 * (maxpwr - tcal[0])) * tcal[1];
	fclose(calfile);
	return(1);
}