Difference between revisions of "File:Pu241heatPlot.jpg"

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define DB double

void ado(FILE *O, int X, int Y)
{ fprintf(O,"%c!PS-Adobe-2.0 EPSF-2.0\n",'%');
        fprintf(O,"%c%cBoundingBox: 0 0 %d %d\n",'%','%',X,Y);
        fprintf(O,"/M {moveto} bind def\n");
        fprintf(O,"/L {lineto} bind def\n");
        fprintf(O,"/S {stroke} bind def\n");
        fprintf(O,"/s {show newpath} bind def\n");
        fprintf(O,"/C {closepath} bind def\n");
        fprintf(O,"/F {fill} bind def\n");
        fprintf(O,"/o {.1 0 360 arc C S} bind def\n");
        fprintf(O,"/times-Roman findfont 20 scalefont setfont\n");
        fprintf(O,"/W {setlinewidth} bind def\n");
        fprintf(O,"/RGB {setrgbcolor} bind def\n");}
//DB casde(DB k, DB p, DB t){ return k*p/(k-p)*(exp(-p*t)-exp(-k*t));}

int main(){int n,m;
DB year=365.25*24*60*60; //second

DB PuT, Puk, Pue;
PuT=14.4 ; // Years
Puk=log(2.)/PuT;
Pue=0.021; //MeV

DB Am, AmT, Amk, Ame;
AmT=432;
Ame=5.638;
Amk=log(2.)/AmT;

DB MeV=1.60218e-13; // Joule
DB Year=365.242374*24*3600; //Second
DB amu=1.660539040e-27 ; //kg

DB k=Puk/Year; //Herz
DB p=Amk/Year; //Herz
DB kp=k*p/(k-p); //Herz
DB W=241.1*amu; //kg , mass of Pu atom
DB E=Ame*MeV; //Joule
DB ko=kp*E; //Watt
DB mo=ko/W; //Watt/kg
printf("kp=%12.4le W=%12.6le E=%12.6le mo=%12.6le\n",kp,W,E,mo);

FILE *o; o=fopen("11.eps","w");
ado(o,1620,1040);
fprintf(o,"10 10 translate 10 10 scale 2 setlinecap 1 setlinejoin\n");
#define M(x,y) fprintf(o,"%6.4lf %6.4lf M\n",0.+x,0.+y);
#define L(x,y) fprintf(o,"%6.4lf %6.4lf L\n",0.+x,0.+y);
M(0,104) L(0,0) L(162,0)
fprintf(o,".5 W S\n");

for(n=10; n<101; n+=10){ M(0,n) L(160,n)}
for(m=10; m<161; m+=10){ M(m,0) L(m,100)}
fprintf(o,".3 W S\n");

DB t;
for(n=0;n<164;n+=2){t=n;
        DB c=exp(-Amk*t)-exp(-Puk*t);
        DB d=c*mo; //Watt/kg
        DB f=d/E; //Bq/kg
        printf("%3.0f %6.4lf %9.3le %9.4lf\n",t,c,f,d);
        if(n==0) M(t,d) else L(t,d)
        }
fprintf(o,"1 0 0 RGB .5 W S\n");

fprintf(o,"showpage\n");
fclose(o);

system("epstopdf 11.eps");
system( "open 11.pdf");

}