home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Club Amiga de Montreal - CAM
/
CAM_CD_1.iso
/
files
/
509.lha
/
ATS_v2.0
/
ode.icn
< prev
next >
Wrap
Text File
|
1991-05-06
|
10KB
|
347 lines
# Ode.icn
# Icon preprocessor to generate c++ code to solve a system of
# ordinary differential equations.
#
# By Dennis J. Darland
# 4/3/91
# Version 1.6
# PROVIDED WITH NO WARRANTY
#
global in,out,n_eq,dep_vars,dep_init_val,x_init_val,h,eqs,iter,here,plot
global lines,ex_eqs,ex_lines,exact,max_terms
procedure main()
dep_vars := []
dep_init_val := []
eqs := []
ex_eqs := []
write("Enter specification file name:")
in_name := read()
write("Enter output c++ file name:")
out_name := read()
in := open(in_name,"r")
out := open(out_name,"w")
write(out,"// ",out_name)
write(out,"// C++ code generated by ode.icn preprocessor for solving")
write(out,"// systems of ordinary differential equations.")
write(out,"// Version 1.6 4/3/91")
write(out,"// PROVIDED WITH NO WARRANTY")
write(out,"#include \"ats.h\"")
line := read(in)
write(out,"// ",line)
line ? (tab(match("exact=")) & exact := tab(many('10')))
line := read(in)
write(out,"// ",line)
line ? (tab(match("max_terms=")) & max_terms := tab(many('1234567890')))
line := read(in)
write(out,"// ",line)
line ? (tab(match("n_eq=")) & n_eq := tab(many('1234567890')))
line := read(in)
write(out,"// ",line)
line ? (tab(match("x=")) & x_init_val := tab(0))
i := 0
while (i<n_eq) do {
line := read(in)
write(out,"// ",line)
line ? (var := tab(upto('=')) & tab(match("=")) & val := tab(0))
dep_vars |||:= [var]
dep_init_val |||:= [val]
i +:= 1
}
line := read(in)
write(out,"// ",line)
line ? (tab(match("h=")) & h := tab(0))
i := 0
while (i<n_eq) do {
lines := []
line := read(in)
write(out,"// ",line)
lines |||:= [line]
while (line ~== "//end") do {
line := read(in)
write(out,"// ",line)
lines |||:= [line]
}
eqs |||:= [lines]
if (exact = 1) then {
ex_lines := []
line := read(in)
write(out,"// ",line)
ex_lines |||:= [line]
while (line ~== "//end") do {
line := read(in)
write(out,"// ",line)
ex_lines |||:= [line]
}
ex_eqs |||:= [ex_lines]
}
i +:= 1
}
line := read(in)
write(out,"// ",line)
line ? (tab(match("iterations=")) & iter := tab(0))
line := read(in)
write(out,"// ",line)
line ? (tab(match("plot=")) & plot := tab(0))
func_decl()
write(out,"// following made global to avoid cfront guru")
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"t_s ",f,", ret_",f,";")
}
write(out,"complex *c_array,rad_xx,p_xx,rad_cp,p_cp,new_h,h_tmp;")
write(out,"double h_re,rad_cp_re,rad_xx_re,p_xx_re,p_cp_re,h_orig_re;")
write(out,"double ex_x,*x_plot,*y_plot[",n_eq,"],*ex_y_plot[",n_eq,"];")
write(out,"static pt_2=0;")
write(out,"FILE *plot_cmds_2;")
write(out,"char fname_2[64];")
write(out,"char cmd_2[81];");
write(out,"main()")
write(out,"\t{")
write(out,"\tregister int i,j;")
write(out,"\tx_plot = (double *) malloc(",iter," * sizeof(double));")
write(out,"\tfor(i=0;i<",n_eq,";i++)")
write(out,"\t\t{");
write(out,"\t\ty_plot[i] = (double *) malloc(",iter," * sizeof(double));")
if (exact = 1) then {
write(out,"\t\tex_y_plot[i] = (double *) malloc(",iter," * sizeof(double));")
}
write(out,"\t\t}");
write(out,"\tc_array= (complex *) malloc(MAX_TERMS * sizeof(complex));")
write(out,"\tfor (i=1;i<MAX_TERMS;i++)")
write(out,"\t\t{")
write(out,"\t\tc_array[i]=complex(0.0,0.0);")
write(out,"\t\t}")
tmp := copy(dep_vars)
tmp_init := copy(dep_init_val)
every f:=gen_pop(tmp) do {
v := pop(tmp_init)
write(out,"\t",f,".set_first(0);")
write(out,"\t",f,".set_last(1);")
write(out,"\tc_array[0]=",v,";")
write(out,"\tset_terms(",f,",c_array);")
write(out,"\t",f,".set_h(",h,");")
write(out,"\th_tmp = ",h,";");
write(out,"\th_orig_re=real(h_tmp);")
write(out,"\t",f,".set_x0(",x_init_val,");")
}
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\t\tret_",f,"=",f,";")
}
write(out,"\tfor (i=0;i<",iter,";i++)")
write(out,"\t\t{")
tmp := copy(dep_vars)
write(out,"\t\tsys_diff_eq(")
f:=gen_pop(tmp)
write(out,"\t\tret_",f,",",f)
every f:=gen_pop(tmp) do {
write(out,"\t\t,ret_",f,",",f)
}
write(out,"\t\t);")
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\t\t",f,"=ret_",f,";")
}
tmp := copy(dep_vars)
if (exact = 1) then {
tmp2 := copy(ex_eqs)
}
write(out,"\t\tj=0;");
every f:=gen_pop(tmp) do {
if (exact = 1) then {
f2 := gen_pop(tmp2)
}
write(out,"\t\tx_plot[i]= get_x0(",f,");");
if (exact = 1) then {
write(out,"\t\tex_x = x_plot[i];")
}
write(out,"\t\tcout << \"",f," is \\n\";")
write(out,"\t\tdisplay(",f,",",plot,");")
write(out,"\t\ty_plot[j][i]=get_y0(",f,");");
if (exact = 1) then {
f3 := gen_pop(f2)
write(out,"\t\tex_y_plot[j][i]=",f3);
every f3 := gen_pop(f2) do {
write(out,f3)
}
write(out,";")
}
write(out,"\t\tj++;");
}
here := 0
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\t\trad_xx = radius_xx(",f,",p_xx);")
write(out,"\t\tcout << \"radius_xx = \" << rad_xx << \" order = \" << p_xx << \" \\n\";")
write(out,"\t\trad_cp = radius_cp(",f,",p_cp);")
write(out,"\t\tcout << \"radius_cp = \" << rad_cp << \" order = \" << p_cp << \" \\n\";")
if (here = 0)
then {
write(out,"\t\tp_xx_re = real(p_xx);")
write(out,"\t\tp_cp_re = real(rad_cp);")
write(out,"\t\trad_xx_re = real(rad_xx);")
write(out,"\t\trad_cp_re = real(rad_cp);")
write(out,"\t\th_re = h_orig_re;")
write(out,"\t\tif (rad_xx_re > 0.0 && p_xx_re > 0.0)")
write(out,"\t\t\t{")
write(out,"\t\t\tif (rad_xx_re / 20.0 < h_re)")
write(out,"\t\t\t\t{")
write(out,"\t\t\t\th_re = rad_xx_re / 20.0;")
write(out,"\t\t\t\t}")
write(out,"\t\t\t}")
write(out,"\t\tif (rad_cp_re > 0.0 && p_cp_re > 0.0)")
write(out,"\t\t\t{")
write(out,"\t\t\tif (rad_cp_re / 20.0 < h_re)")
write(out,"\t\t\t\t{")
write(out,"\t\t\t\th_re = rad_cp_re / 20.0;")
write(out,"\t\t\t\t}")
write(out,"\t\t\t}")
write(out,"\t\tnew_h = complex(h_re,0.0);")
}
write(out,"\t\tadjust_h(",f,",new_h);")
here := 1
}
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\t\tjump(",f,");")
}
write(out,"\t\t}")
write(out,"\tfor(pt_2 = 0;pt_2 <",n_eq,";pt_2++)")
write(out,"\t\t{")
write(out,"\t\tsprintf(fname_2,\"ram:mapleplot2_%d\",pt_2);")
write(out,"\t\tplot_cmds_2 = fopen(\"ram:plottmp\",\"w\");")
write(out,"\t\tfprintf(plot_cmds_2,\"plotdevice := amiga;\\n\");")
write(out,"\t\tfprintf(plot_cmds_2,\"plot([\");")
write(out,"\t\tfor(i=0;i<",iter,";i++)")
write(out,"\t\t\t{");
write(out,"\t\t\tif (i != 0)")
write(out,"\t\t\t\tfprintf(plot_cmds_2,\",\");")
write(out,"\t\t\tfprintf(plot_cmds_2,\"%g,%g\\n\",")
write(out,"\t\t\tx_plot[i],y_plot[pt_2][i]);")
write(out,"\t\t\t}")
write(out,"\t\tfprintf(plot_cmds_2,\"]);\\n\");")
write(out,"\t\tfprintf(plot_cmds_2,\"quit\\n\");")
write(out,"\t\tfclose(plot_cmds_2);")
write(out,"\t\tsprintf(cmd_2,\"iconz <ram:plottmp >%s pltcvt\",fname_2);")
write(out,"\t\tsystem(cmd_2);")
if (plot ~= 0) then {
write(out,"\t\tsprintf(cmd_2,\"maple <%s\",fname_2);")
write(out,"\t\tsystem(cmd_2);")
}
if (exact = 1) then {
write(out,"\t\tsprintf(fname_2,\"ram:mapleplot_ex_%d\",pt_2);")
write(out,"\t\tplot_cmds_2 = fopen(\"ram:plottmp\",\"w\");")
write(out,"\t\tfprintf(plot_cmds_2,\"plotdevice := amiga;\\n\");")
write(out,"\t\tfprintf(plot_cmds_2,\"plot([\");")
write(out,"\t\tfor(i=0;i<",iter,";i++)")
write(out,"\t\t\t{");
write(out,"\t\t\tif (i != 0)")
write(out,"\t\t\t\tfprintf(plot_cmds_2,\",\");")
write(out,"\t\t\tfprintf(plot_cmds_2,\"%g,%g\\n\",")
write(out,"\t\t\tx_plot[i],ex_y_plot[pt_2][i]);")
write(out,"\t\t\t}")
write(out,"\t\tfprintf(plot_cmds_2,\"]);\\n\");")
write(out,"\t\tfprintf(plot_cmds_2,\"quit\\n\");")
write(out,"\t\tfclose(plot_cmds_2);")
write(out,"\t\tsprintf(cmd_2,\"iconz <ram:plottmp >%s pltcvt\",fname_2);")
write(out,"\t\tsystem(cmd_2);")
if (plot ~= 0) then {
write(out,"\t\tsprintf(cmd_2,\"maple <%s\",fname_2);")
write(out,"\t\tsystem(cmd_2);")
}
}
write(out,"\t\t}")
write(out,"\t}")
# end of main()
write(out,"void equation (t_s &dy,t_s &y,...)")
write(out,"\t{ // function to resolve extern in ats.cp")
write(out,"\t}")
tmp := copy(dep_vars)
tmp_eqs := copy(eqs)
every f:=gen_pop(tmp) do {
eq := pop(tmp_eqs)
eq2 := copy(eq)
tmp2 := copy(dep_vars)
write(out,"void equation_",f,"(t_s &diff_",f,)
every f2:=gen_pop(tmp2) do {
write(out,",t_s &",f2)
}
write(out,")")
write(out,"\t{")
every l := gen_pop(eq2) do {
write(out,"\t",l)
}
write(out,"\t;")
write(out,"\t}")
}
tmp := copy(dep_vars)
write(out,"void sys_diff_eq(")
f:=gen_pop(tmp)
write(out,"t_s &ret_",f,",t_s &",f)
every f:=gen_pop(tmp) do {
write(out,",t_s &ret_",f,",t_s &",f)
}
write(out,")")
tmp:=copy(dep_vars)
f := pop(tmp)
write(out,"\t{")
write(out,"\tregister int k;")
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\tt_s diff_",f,";")
}
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\tdiff_",f,"=",f,";")
}
write(out,"\tfor (k = get_last_term(",f,");k<",max_terms,";k++)")
write(out,"\t\t{")
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
tmp2 := copy(dep_vars)
write(out,"\t\tequation_",f,"(diff_",f,)
every f2:=gen_pop(tmp2) do {
write(out,"\t\t,",f2)
}
write(out,"\t\t);")
}
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\t\tset_term_from_deriv(k,diff_",f,",",f,");")
}
write(out,"\t\t}")
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
write(out,"\tret_",f,"=",f,";")
}
write(out,"\t}")
end
procedure func_decl()
tmp := copy(dep_vars)
write(out,"void sys_diff_eq(")
f:=gen_pop(tmp)
write(out,"t_s &ret_",f,",t_s &",f)
every f:=gen_pop(tmp) do {
write(out,",t_s &ret_",f,",t_s &",f)
}
write(out,");")
tmp := copy(dep_vars)
every f:=gen_pop(tmp) do {
tmp2 := copy(dep_vars)
write(out,"void equation_",f,"(t_s &diff_",f,)
every f2:=gen_pop(tmp2) do {
write(out,",t_s &",f2)
}
write(out,");")
}
end
procedure gen_pop(l)
if (t := pop(l))
then suspend t | gen_pop(l)
else fail
end