diff -urN ORG/lib/@mkall lib/@mkall
--- ORG/lib/@mkall	2017-12-23 18:00:00.000000000 +0900
+++ lib/@mkall	2021-02-22 12:22:06.000000000 +0900
@@ -4,6 +4,7 @@
 #            to generate Object xxxxx.o
 #   'fort' must be aliased to the desired FORTRAN compiler
 #---------------------------------------
+gcc  -Wall -c axis.c
 gcc  -Wall -c calma.c
 gcc  -Wall -c cont.c
 fort -Wall -c getargs.f90
diff -urN ORG/lib/axis.c lib/axis.c
--- ORG/lib/axis.c	1970-01-01 09:00:00.000000000 +0900
+++ lib/axis.c	2021-02-22 12:22:06.000000000 +0900
@@ -0,0 +1,199 @@
+#include    <stdio.h>
+#include    <stdlib.h>
+#include    <math.h>
+
+int plot(double xp, double yp, int md);
+int where(double *wx, double *wy, double *wfac);
+int pcstr(double xp, double yp, double chi, char *text, double ang, int ns);
+int pmark(int mark, double xc, double yc, double size, double thick, int mcol);
+
+#define RAD     (180./3.14159)
+
+/*************/
+/*   lmark   */
+/*************/
+
+static int mks[10] = { 2, 0, 4, 8, 9, 3, 1, 5, 6, 7 };
+int lmark(int lmk, double size)
+  // ƒyƒ“‚ÌŒ»�ÝˆÊ’u‚Épmarkƒ‹�[ƒ`ƒ“‚Æ“¯—l‚Ìƒ}�[ƒN‚ð�•�ü‚Å•`‚­‚ª�C
+  // ‘æ1ˆø�” lmk ‚Ì’l‚Ì 0�`9 ‚É‘Î‰ž‚µ‚Ä�Cpmarkƒ‹�[ƒ`ƒ“‚Ì‘æ1ˆø�” mark ‚ª
+  // 2, 0, 4, 8, 9, 3, 1, 5, 6, 7 ‚Ìƒ}�[ƒN‚ð•`‚­�D
+  // ’A‚µ�Clmk ‚É •‰ ‚Ü‚½‚Í 10 ˆÈ�ã‚Ì’l‚ð—^‚¦‚é‚Æ�C�¬‚³‚¢ (Žw’èsize ‚Ì 1/3)
+  // ‚Ê‚è‚Â‚Ô‚µ‚½�•ŠÛ‚ð•`‚­�D
+{   double wx, wy, wfac;
+    if (where(&wx,&wy,&wfac) < 0)
+        { fprintf(stderr, "lmark: error at calling where\n"); return(-1); }
+    if ((lmk<0) || (lmk>9)) {
+        return(pmark(0, wx, wy, size/3., 0., 0));
+    } else {
+        return(pmark(mks[lmk], wx, wy, size, size*0.075, 0));
+    }
+}
+
+/***********************/
+/*   symbol / number   */
+/***********************/
+
+int symbol(double xp, double yp, double high, char *text, double ang, int nc)
+{   double w, wf;
+    int k;
+    if (xp == 999.) { where(&xp, &w, &wf); if (yp == 999.) yp = w; }
+    else if (yp == 999.) where(&w, &yp, &wf);
+    if (nc <= 0) {
+        if (nc < -1) k = plot(xp,yp,2); else k = plot(xp,yp,3);
+        if (k < 0) return(-1);
+        lmark((int)*text, high);
+        //    [line-color and thickness are fixed to black and high*0.075]
+    }
+    else {
+        if (plot(xp,yp,3) < 0) return(-1);
+        pcstr(xp, yp, high, text, ang, nc);
+    }
+    return(0);
+}
+
+int number(double xp, double yp, double high, double fv, double ang, int ndec)
+{   double w, wf;
+    char str[32], fmt[8], *s;
+    int ns;
+    if (ndec < -9) ndec = -9;
+    if (ndec > +9) ndec = +9;
+    if (ndec > 0) sprintf(fmt, "%%%2d.%df", 22+ndec, ndec);
+    else if (ndec == 0) sprintf(fmt, "%%21.0f.");
+    else { while (++ndec < 0) fv /= 10.; sprintf(fmt, "%%21.0f"); }
+    if (sprintf(str,fmt,fv) < 0)
+        { fprintf(stderr, "number: too long expression\n"); return(2); }
+    for (s=str; *s==' '; s++) if (*s == '\0') return(1);
+    if (xp == 999.) { where(&xp, &w, &wf); if (yp == 999.) yp = w; }
+    else if (yp == 999.) where(&w, &yp, &wf);
+    if (plot(xp,yp,3) < 0) return(-1);
+    ns = &str[31]+1 - s;  pcstr(xp, yp, high, s, ang, ns);
+    return(0);
+}
+
+/********************/
+/*   scale / line   */
+/********************/
+
+int scale(float va[], double axlen, int npt, int inc)
+{   int n, ipt, ke, kf;
+    double vl, vh, vd, ul, uh, uf;
+    if ((npt<=1) || (inc<=0) || (axlen<2.))
+        { fprintf(stderr, "scale: parameter value invalid\n"); return(-1); }
+    vl = vh = va[0];
+    for (n=inc,ipt=1; ipt<npt; ipt++,n+=inc)
+        { if (va[n] < vl) vl = va[n]; if (va[n] > vh) vh = va[n]; }
+    if (vl == vh)
+        { fprintf(stderr, "scale: identical data values\n"); return(-1); }
+    vd = (vh-vl)/axlen; ke = (int)log10(vd); if (vd < 1.) ke--;
+    kf = (int)(vd/pow(10.,(double)ke));
+    do {
+        kf++;
+        if (kf > 8) { kf = 1; ke++; }
+        else if (kf > 5) kf = 8;
+        else if (kf == 3) kf = 4;
+        vd = (double)kf*pow(10.,(double)ke);
+        ul = vl/vd; uh = vh/vd;
+        uf = (double)(int)ul; if (ul < uf) uf = uf-1.;
+    } while ((uh-uf) > axlen);
+    va[n] = uf*vd; va[n+inc] = vd;
+    return(0);
+}
+
+int line(float xa[], float ya[], int npt, int inc, int ltype, int lmk)
+{   int nf, n, ipt;
+    double xf, yf, xd, yd, x, y;
+    if ((npt<=0) || (inc<=0))
+        { fprintf(stderr, "line: parameter value invalid\n"); return(-1); }
+    nf = npt*inc; xf = xa[nf]; yf = ya[nf]; xd = xa[nf+inc]; yd = ya[nf+inc];
+    if (xd == 0.) xd = 1.;
+    if (yd == 0.) yd = 1.;
+    x = (xa[0]-xf)/xd; y = (ya[0]-yf)/yd; if (plot(x,y,3) < 0) return(-1);
+    if (ltype != 0) lmark(lmk, 0.21);
+    for (n=inc,ipt=1; ipt<npt; ipt++,n+=inc) {
+        x = (xa[n]-xf)/xd; y = (ya[n]-yf)/yd;
+        if (ltype < 0) plot(x, y, 3); else plot(x, y, 2);
+        if (ltype != 0) {
+            if ((ipt%ltype) == 0) lmark(lmk, 0.21);
+        }
+    }
+    return(0);
+}
+
+/************/
+/*   axis   */
+/************/
+
+int axis(double xp, double yp, char *text, int ltx,
+         double axlen, double angle, double firstv, double deltav)
+{   int ke, i;
+    double bt, bv, bs, ca, sa, fe, as, xs, ys, xt, yt, vf;
+    if ((ltx==0) || (axlen<2.))
+        { fprintf(stderr, "axis: parameter value invalid\n"); return(-1); }
+    if (ltx > 0)     { bt =  0.735; bv =  0.298; bs =  0.178; }
+    else { ltx = -ltx; bt = -1.015; bv = -0.508; bs = -0.178; }
+    ca = cos(angle/RAD); sa = sin(angle/RAD);
+    if (deltav == 0.) deltav = 1.;
+    ke = (int)log10(fabs(deltav*1.00001));
+    if (ke > 0) ke--;
+    else if (ke < 0) ke++;
+    fe = pow(10.,(double)ke); firstv /= fe; deltav /= fe;
+    as = axlen/2.-0.14*(double)ltx; if (ke != 0) as -= 0.88;
+    xs = xp+ca*as-sa*bt; ys = yp+sa*as+ca*bt;
+    if (plot(xs,ys,3) < 0) return(-1);
+    pcstr(xs, ys, 0.28, text, angle, ltx);
+    if (ke != 0) {
+        as = ltx*0.28; xs += ca*as; ys += sa*as;
+        pcstr(xs, ys, 0.28, "  *10", angle, 5);
+        as = 5*0.28; xs += ca*as-sa*0.18; ys += sa*as+ca*0.18;
+        number(xs, ys, 0.18, (double)ke, angle, -1);
+    }
+    xt = xp; yt = yp; plot(xt, yt, 3);
+    xs = xt-sa*bs; ys = yt+ca*bs; plot(xs, ys, 2);
+    xs = xt-ca*0.28-sa*bv; ys = yt-sa*0.28+ca*bv; vf = firstv;
+    number(xs, ys, 0.21, vf, angle, 2);
+    for (i=1; i<=(int)axlen; i++) {
+        xt += ca; yt += sa; vf += deltav; plot(xt, yt, 3);
+        xs = xt-sa*bs; ys = yt+ca*bs; plot(xs, ys, 2);
+        if ((i%2) == 0) {
+            xs = xt-ca*0.49-sa*bv; ys = yt-sa*0.49+ca*bv;
+            number(xs, ys, 0.21, vf, angle, 2);
+        }
+    }
+    xs = xp+ca*axlen; ys = yp+sa*axlen; plot(xs, ys, 3); plot(xp, yp, 2);
+    return(0);
+}
+
+/*******************************************/
+/*   FORTRAN Language Interface Routines   */
+/*******************************************/
+
+#define W       double
+
+/***********************************************************/
+/*   lmark_ / symbol_ / number_ / scale_ / line_ / axis_   */
+/***********************************************************/
+
+
+void lmark_(int *lmk, float *size)
+{ if (lmark(*lmk,(W)*size) < 0) exit(1); }
+
+void symbol_(float *xp, float *yp, float *high, char *text,
+             float *ang, int *nc)
+{ if (symbol((W)*xp,(W)*yp,(W)*high,text,(W)*ang,*nc) < 0) exit(1); }
+
+void number_(float *xp, float *yp, float *high, float *fv,
+             float *ang, int *ndec)
+{ if (number((W)*xp,(W)*yp,(W)*high,(W)*fv,(W)*ang,*ndec) < 0) exit(1); }
+
+void scale_(float va[], float *axlen, int *npt, int *inc)
+{ if (scale(va,(W)*axlen,*npt,*inc) < 0) exit(1); }
+
+void line_(float xa[], float ya[], int *npt, int *inc, int *ltype, int *mark)
+{ if (line(xa,ya,*npt,*inc,*ltype,*mark) < 0) exit(1); }
+
+void axis_(float *xp, float *yp, char *text, int *ltx,
+           float *axlen, float *angle, float *firstv, float *deltav)
+{ if (axis((W)*xp,(W)*yp,text,*ltx,(W)*axlen,
+           (W)*angle,(W)*firstv,(W)*deltav) < 0) exit(1); }
+
