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 >
Encoding:
Text File  |  1991-05-07  |  9.6 KB  |  347 lines
  1. # Ode.icn
  2. # Icon preprocessor to generate c++ code to solve a system of
  3. # ordinary differential equations. 
  4. #
  5. # By Dennis J. Darland
  6. # 4/3/91
  7. # Version 1.6
  8. # PROVIDED WITH NO WARRANTY
  9. #
  10. global in,out,n_eq,dep_vars,dep_init_val,x_init_val,h,eqs,iter,here,plot
  11. global lines,ex_eqs,ex_lines,exact,max_terms
  12. procedure main()
  13. dep_vars := []
  14. dep_init_val := []
  15. eqs := []
  16. ex_eqs := []
  17. write("Enter specification file name:")
  18. in_name := read()
  19. write("Enter output c++ file name:")
  20. out_name := read()
  21. in := open(in_name,"r")
  22. out := open(out_name,"w")
  23. write(out,"// ",out_name)
  24. write(out,"// C++ code generated by ode.icn preprocessor for solving")
  25. write(out,"// systems of ordinary differential equations.")
  26. write(out,"// Version 1.6 4/3/91")
  27. write(out,"// PROVIDED WITH NO WARRANTY")
  28. write(out,"#include \"ats.h\"")
  29. line := read(in)
  30. write(out,"// ",line)
  31. line ? (tab(match("exact=")) & exact := tab(many('10')))
  32. line := read(in)
  33. write(out,"// ",line)
  34. line ? (tab(match("max_terms=")) & max_terms := tab(many('1234567890')))
  35. line := read(in)
  36. write(out,"// ",line)
  37. line ? (tab(match("n_eq=")) & n_eq := tab(many('1234567890')))
  38. line := read(in)
  39. write(out,"// ",line)
  40. line ? (tab(match("x=")) & x_init_val := tab(0))
  41. i := 0
  42. while (i<n_eq) do {
  43.     line := read(in)
  44.     write(out,"// ",line)
  45.     line ? (var := tab(upto('=')) & tab(match("=")) & val := tab(0))
  46.     dep_vars |||:= [var]
  47.     dep_init_val |||:= [val]
  48.     i +:= 1
  49.     }
  50. line := read(in)
  51. write(out,"// ",line)
  52. line ? (tab(match("h=")) & h := tab(0))
  53. i := 0
  54. while (i<n_eq) do {
  55.     lines := []
  56.     line := read(in)
  57.     write(out,"// ",line)
  58.     lines |||:= [line]
  59.     while (line ~== "//end") do {
  60.         line := read(in)
  61.         write(out,"// ",line)
  62.         lines |||:= [line]
  63.         }
  64.     eqs |||:= [lines]
  65.     if (exact = 1) then {
  66.         ex_lines := []
  67.         line := read(in)
  68.         write(out,"// ",line)
  69.         ex_lines |||:= [line]
  70.         while (line ~== "//end") do {
  71.             line := read(in)
  72.             write(out,"// ",line)
  73.             ex_lines |||:= [line]
  74.             }
  75.         ex_eqs |||:= [ex_lines]
  76.         }
  77.     i +:= 1
  78.     }
  79. line := read(in)
  80. write(out,"// ",line)
  81. line ? (tab(match("iterations=")) & iter := tab(0))
  82. line := read(in)
  83. write(out,"// ",line)
  84. line ? (tab(match("plot=")) & plot := tab(0))
  85. func_decl()
  86. write(out,"// following made global to avoid cfront guru")
  87. tmp := copy(dep_vars)
  88.  every f:=gen_pop(tmp) do {
  89.      write(out,"t_s ",f,", ret_",f,";")
  90.     }
  91. write(out,"complex *c_array,rad_xx,p_xx,rad_cp,p_cp,new_h,h_tmp;")
  92. write(out,"double h_re,rad_cp_re,rad_xx_re,p_xx_re,p_cp_re,h_orig_re;")
  93. write(out,"double ex_x,*x_plot,*y_plot[",n_eq,"],*ex_y_plot[",n_eq,"];")
  94. write(out,"static pt_2=0;")
  95. write(out,"FILE *plot_cmds_2;")
  96. write(out,"char fname_2[64];")
  97. write(out,"char cmd_2[81];");
  98. write(out,"main()")
  99. write(out,"\t{")
  100. write(out,"\tregister int i,j;")
  101. write(out,"\tx_plot = (double *) malloc(",iter," * sizeof(double));")
  102. write(out,"\tfor(i=0;i<",n_eq,";i++)")
  103. write(out,"\t\t{");
  104. write(out,"\t\ty_plot[i] = (double *) malloc(",iter," * sizeof(double));")
  105. if (exact = 1) then {
  106. write(out,"\t\tex_y_plot[i] = (double *) malloc(",iter," * sizeof(double));")
  107.     }
  108. write(out,"\t\t}");
  109. write(out,"\tc_array= (complex *) malloc(MAX_TERMS * sizeof(complex));")
  110. write(out,"\tfor (i=1;i<MAX_TERMS;i++)")
  111. write(out,"\t\t{")
  112. write(out,"\t\tc_array[i]=complex(0.0,0.0);")
  113. write(out,"\t\t}")
  114. tmp := copy(dep_vars)
  115. tmp_init := copy(dep_init_val)
  116.  every f:=gen_pop(tmp) do {
  117.      v := pop(tmp_init)
  118.      write(out,"\t",f,".set_first(0);")
  119.      write(out,"\t",f,".set_last(1);")
  120.     write(out,"\tc_array[0]=",v,";")
  121.      write(out,"\tset_terms(",f,",c_array);")
  122.     write(out,"\t",f,".set_h(",h,");")
  123.     write(out,"\th_tmp = ",h,";");
  124.     write(out,"\th_orig_re=real(h_tmp);")
  125.     write(out,"\t",f,".set_x0(",x_init_val,");")
  126.     }
  127. tmp := copy(dep_vars)
  128.  every f:=gen_pop(tmp) do {
  129.      write(out,"\t\tret_",f,"=",f,";")
  130.     }
  131. write(out,"\tfor (i=0;i<",iter,";i++)")
  132. write(out,"\t\t{")
  133. tmp := copy(dep_vars)
  134. write(out,"\t\tsys_diff_eq(")
  135.  f:=gen_pop(tmp)
  136.  write(out,"\t\tret_",f,",",f)
  137.  every f:=gen_pop(tmp) do {
  138.      write(out,"\t\t,ret_",f,",",f)
  139.     }
  140. write(out,"\t\t);")
  141. tmp := copy(dep_vars)
  142.  every f:=gen_pop(tmp) do {
  143.      write(out,"\t\t",f,"=ret_",f,";")
  144.     }
  145. tmp := copy(dep_vars)
  146. if (exact = 1) then {
  147.     tmp2 := copy(ex_eqs)
  148.     }
  149. write(out,"\t\tj=0;");
  150.  every f:=gen_pop(tmp) do {
  151.      if (exact = 1) then {
  152.     f2 := gen_pop(tmp2)
  153.         }
  154.     write(out,"\t\tx_plot[i]= get_x0(",f,");");
  155.     if (exact = 1) then {
  156.         write(out,"\t\tex_x = x_plot[i];")
  157.         }
  158.      write(out,"\t\tcout << \"",f," is \\n\";") 
  159.      write(out,"\t\tdisplay(",f,",",plot,");")
  160.     write(out,"\t\ty_plot[j][i]=get_y0(",f,");");
  161.     if (exact = 1) then {
  162.         f3 := gen_pop(f2)
  163.         write(out,"\t\tex_y_plot[j][i]=",f3);
  164.         every f3 := gen_pop(f2) do {
  165.              write(out,f3)
  166.             }
  167.         write(out,";")
  168.         } 
  169.     write(out,"\t\tj++;");
  170.     }
  171. here := 0
  172. tmp := copy(dep_vars)
  173.  every f:=gen_pop(tmp) do {
  174.      write(out,"\t\trad_xx = radius_xx(",f,",p_xx);")
  175.     write(out,"\t\tcout << \"radius_xx = \" << rad_xx << \" order = \" << p_xx << \" \\n\";") 
  176.      write(out,"\t\trad_cp = radius_cp(",f,",p_cp);")
  177.     write(out,"\t\tcout << \"radius_cp = \" << rad_cp << \" order = \" << p_cp << \" \\n\";") 
  178.     if (here = 0) 
  179.     then {
  180.         write(out,"\t\tp_xx_re = real(p_xx);")
  181.         write(out,"\t\tp_cp_re = real(rad_cp);")
  182.         write(out,"\t\trad_xx_re = real(rad_xx);")
  183.         write(out,"\t\trad_cp_re = real(rad_cp);")
  184.         write(out,"\t\th_re = h_orig_re;")
  185.         write(out,"\t\tif (rad_xx_re > 0.0 && p_xx_re > 0.0)")
  186.         write(out,"\t\t\t{")
  187.         write(out,"\t\t\tif (rad_xx_re / 20.0 < h_re)")
  188.         write(out,"\t\t\t\t{")
  189.         write(out,"\t\t\t\th_re = rad_xx_re / 20.0;")
  190.         write(out,"\t\t\t\t}")
  191.         write(out,"\t\t\t}")
  192.         write(out,"\t\tif (rad_cp_re > 0.0 && p_cp_re > 0.0)")
  193.         write(out,"\t\t\t{")
  194.         write(out,"\t\t\tif (rad_cp_re / 20.0 < h_re)")
  195.         write(out,"\t\t\t\t{")
  196.         write(out,"\t\t\t\th_re = rad_cp_re / 20.0;")
  197.         write(out,"\t\t\t\t}")
  198.         write(out,"\t\t\t}")
  199.         write(out,"\t\tnew_h = complex(h_re,0.0);") 
  200.         }
  201.     write(out,"\t\tadjust_h(",f,",new_h);")
  202.      here := 1
  203.     }
  204.  
  205. tmp := copy(dep_vars)
  206.  every f:=gen_pop(tmp) do {
  207.      write(out,"\t\tjump(",f,");")
  208.     }
  209. write(out,"\t\t}")
  210. write(out,"\tfor(pt_2 = 0;pt_2 <",n_eq,";pt_2++)") 
  211. write(out,"\t\t{")
  212. write(out,"\t\tsprintf(fname_2,\"ram:mapleplot2_%d\",pt_2);")
  213. write(out,"\t\tplot_cmds_2 = fopen(\"ram:plottmp\",\"w\");")
  214. write(out,"\t\tfprintf(plot_cmds_2,\"plotdevice := amiga;\\n\");")
  215. write(out,"\t\tfprintf(plot_cmds_2,\"plot([\");")
  216. write(out,"\t\tfor(i=0;i<",iter,";i++)")
  217. write(out,"\t\t\t{");
  218. write(out,"\t\t\tif (i != 0)")
  219. write(out,"\t\t\t\tfprintf(plot_cmds_2,\",\");")
  220. write(out,"\t\t\tfprintf(plot_cmds_2,\"%g,%g\\n\",")
  221. write(out,"\t\t\tx_plot[i],y_plot[pt_2][i]);")
  222. write(out,"\t\t\t}")
  223. write(out,"\t\tfprintf(plot_cmds_2,\"]);\\n\");")
  224. write(out,"\t\tfprintf(plot_cmds_2,\"quit\\n\");")
  225. write(out,"\t\tfclose(plot_cmds_2);")
  226. write(out,"\t\tsprintf(cmd_2,\"iconz <ram:plottmp >%s pltcvt\",fname_2);")
  227. write(out,"\t\tsystem(cmd_2);")
  228. if (plot ~= 0) then {
  229.     write(out,"\t\tsprintf(cmd_2,\"maple <%s\",fname_2);")
  230.     write(out,"\t\tsystem(cmd_2);")
  231.     }
  232. if (exact = 1) then {
  233. write(out,"\t\tsprintf(fname_2,\"ram:mapleplot_ex_%d\",pt_2);")
  234. write(out,"\t\tplot_cmds_2 = fopen(\"ram:plottmp\",\"w\");")
  235. write(out,"\t\tfprintf(plot_cmds_2,\"plotdevice := amiga;\\n\");")
  236. write(out,"\t\tfprintf(plot_cmds_2,\"plot([\");")
  237. write(out,"\t\tfor(i=0;i<",iter,";i++)")
  238. write(out,"\t\t\t{");
  239. write(out,"\t\t\tif (i != 0)")
  240. write(out,"\t\t\t\tfprintf(plot_cmds_2,\",\");")
  241. write(out,"\t\t\tfprintf(plot_cmds_2,\"%g,%g\\n\",")
  242. write(out,"\t\t\tx_plot[i],ex_y_plot[pt_2][i]);")
  243. write(out,"\t\t\t}")
  244. write(out,"\t\tfprintf(plot_cmds_2,\"]);\\n\");")
  245. write(out,"\t\tfprintf(plot_cmds_2,\"quit\\n\");")
  246. write(out,"\t\tfclose(plot_cmds_2);")
  247. write(out,"\t\tsprintf(cmd_2,\"iconz <ram:plottmp >%s pltcvt\",fname_2);")
  248. write(out,"\t\tsystem(cmd_2);")
  249. if (plot ~= 0) then {
  250.     write(out,"\t\tsprintf(cmd_2,\"maple <%s\",fname_2);")
  251.     write(out,"\t\tsystem(cmd_2);")
  252.     }
  253.     }
  254. write(out,"\t\t}")
  255. write(out,"\t}")
  256. # end of main()
  257. write(out,"void equation (t_s &dy,t_s &y,...)")
  258. write(out,"\t{ // function to resolve extern in ats.cp")
  259. write(out,"\t}")
  260. tmp := copy(dep_vars)
  261. tmp_eqs := copy(eqs)
  262.  every f:=gen_pop(tmp) do {
  263.      eq := pop(tmp_eqs)
  264.     eq2 := copy(eq) 
  265.      tmp2 := copy(dep_vars)
  266.     write(out,"void equation_",f,"(t_s &diff_",f,)
  267.      every f2:=gen_pop(tmp2) do {
  268.      write(out,",t_s &",f2)
  269.      }
  270.     write(out,")")
  271.     write(out,"\t{")
  272.     every l := gen_pop(eq2) do {
  273.         write(out,"\t",l)
  274.         }
  275.     write(out,"\t;")
  276.     write(out,"\t}")
  277.     }
  278. tmp := copy(dep_vars)
  279. write(out,"void sys_diff_eq(")
  280.    f:=gen_pop(tmp)
  281.  write(out,"t_s &ret_",f,",t_s &",f)
  282.  every f:=gen_pop(tmp) do {
  283.      write(out,",t_s &ret_",f,",t_s &",f)
  284.     }
  285. write(out,")")
  286. tmp:=copy(dep_vars)
  287. f := pop(tmp)
  288. write(out,"\t{")
  289. write(out,"\tregister int k;")
  290. tmp := copy(dep_vars)
  291.  every f:=gen_pop(tmp) do {
  292.      write(out,"\tt_s diff_",f,";")
  293.     }
  294. tmp := copy(dep_vars)
  295.  every f:=gen_pop(tmp) do {
  296.      write(out,"\tdiff_",f,"=",f,";")
  297.     }
  298. write(out,"\tfor (k = get_last_term(",f,");k<",max_terms,";k++)")
  299. write(out,"\t\t{")
  300. tmp := copy(dep_vars)
  301.  every f:=gen_pop(tmp) do {
  302.      tmp2 := copy(dep_vars)
  303.     write(out,"\t\tequation_",f,"(diff_",f,)
  304.      every f2:=gen_pop(tmp2) do {
  305.      write(out,"\t\t,",f2)
  306.      }
  307.     write(out,"\t\t);")
  308.     }
  309. tmp := copy(dep_vars)
  310.  every f:=gen_pop(tmp) do {
  311.      write(out,"\t\tset_term_from_deriv(k,diff_",f,",",f,");")
  312.     }
  313. write(out,"\t\t}")
  314. tmp := copy(dep_vars)
  315.  every f:=gen_pop(tmp) do {
  316.      write(out,"\tret_",f,"=",f,";")
  317.     }
  318. write(out,"\t}")
  319. end
  320.  
  321. procedure func_decl()
  322. tmp := copy(dep_vars)
  323. write(out,"void sys_diff_eq(")
  324.    f:=gen_pop(tmp)
  325.  write(out,"t_s &ret_",f,",t_s &",f)
  326.  every f:=gen_pop(tmp) do {
  327.      write(out,",t_s &ret_",f,",t_s &",f)
  328.     }
  329. write(out,");")
  330.  
  331. tmp := copy(dep_vars)
  332.  every f:=gen_pop(tmp) do {
  333.      tmp2 := copy(dep_vars)
  334.     write(out,"void equation_",f,"(t_s &diff_",f,)
  335.      every f2:=gen_pop(tmp2) do {
  336.      write(out,",t_s &",f2)
  337.      }
  338.     write(out,");")
  339.     }
  340. end
  341.  
  342. procedure gen_pop(l)
  343.     if (t := pop(l))
  344.      then suspend t | gen_pop(l)
  345.      else fail
  346. end        
  347.