program slope;
uses crt,graph;
const
 nx=100;mx=4;my=4;
var
 maxx,maxy,graphdriver,graphmode,errorcode:integer;
 xasp,yasp:word;xa,ya,xb,yb,dx,dy,sx,sy:real;
 l,h,h1,h2,gw,gd1,gn1,cc1,phi1,gd2,gn2,cc2,phi2,x1,y1,x2,y2,y0:real;
 s,p:array[1..nx] of real;ff:array[0..mx,0..my] of real;data:text;
procedure title;
 begin
  clrscr;gotoxy(36,1);textbackground(7);textcolor(0);write(' SLOPE ');
  textbackground(0);textcolor(7);writeln;writeln;
 end;
procedure graphinitialize;
 begin
  graphdriver:=detect;initgraph(graphdriver,graphmode,'');
  errorcode:=graphresult;
  if (errorcode<>grok) then
   begin
    writeln('Error in graphics :',grapherrormsg(errorcode));
    writeln;writeln('Program interrupted.');halt(1);
   end;
  setcolor(7);setbkcolor(0);setlinestyle(0,0,1);
  setfillstyle(11,7);maxx:=getmaxx;maxy:=getmaxy;
  getaspectratio(xasp,yasp);closegraph;
 end;
procedure dots(x1,y1,x2,y2:integer);
 var
  i,x,y,nr:integer;dx,dy,dr,xa,ya:real;
 begin
  xa:=x1;dx:=x2-x1;ya:=y1;dy:=y2-y1;
  if dx<0 then begin xa:=x2;ya:=y2;dx:=-dx;dy:=-dy;end;
  if (x1=x2) and (dy<0) then
   begin
    xa:=x2;ya:=y2;dx:=-dx;dy:=-dy;
   end;
  dr:=(sqrt(dx*dx+4*dy*dy))/4;
  if (dr>0) then begin dx:=dx/dr;dy:=dy/dr;end;
  nr:=trunc(dr);
  for i:=0 to nr do
   begin
    x:=round(xa);y:=round(ya);
    line(x,y,x,y);xa:=xa+dx;ya:=ya+dy;
   end;
 end;
procedure input;
 var
  name:string;
 begin
  title;
  writeln('This is a program for the analysis of the stability');
  writeln('of a slope, using the Bishop method.');writeln;
  write('Name of input datafile ............ ');readln(name);
  assign(data,name);reset(data);readln(data,l,h,h1,h2,gw);
  readln(data,gd1,gn1,cc1,phi1);readln(data,gd2,gn2,cc2,phi2);
  readln(data,x1,y1,x2,y2,y0);
  if (y1<h) then y1:=h;if (y2<y1+0.005) then y2:=y1+0.005;
  if (x2<x1+0.005) then x2:=x1+0.005;
  if (h2>h) then h2:=h;if (h1>h) then h1:=h;
  if (h1<0.0) then h1:=0.0;if (h2<0.0) then h1:=h2;
  close(data);title;
 end;
procedure stability;
 var
  i,j,k,kk,ia,ib,ja,jb:integer;
  xc,yc,xl,yl,xr,yr,r,f,fa,pi,cc,ph,tf,a,b,bb:real;
  dx,x,x3,x4,y,yb,yt,yn,co,si,ta,tb,e,p1,p2,p3:real;
 begin
  clrscr;
  xa:=-l;xb:=2*l;ya:=-h;yb:=h;
  if (x1<xa) then xa:=x1;if (x1>xb) then xb:=x1;
  if (y1<ya) then ya:=y1;if (y1>yb) then yb:=y1;
  if (x2<xa) then xa:=x2;if (x2>xb) then xb:=x2;
  if (y2<ya) then ya:=y2;if (y2>yb) then yb:=y2;
  if (y0<ya) then ya:=y0;if (y0>yb) then yb:=y0;
  dx:=xb-xa;dy:=yb-ya;sx:=maxx/dx;sy:=(yasp/xasp)*maxy/dy;
  if sy<sx then sx:=sy;sy:=xasp*sx/yasp;
  pi:=3.1415926;
  for i:=0 to mx do
   begin
    xc:=x1+i*(x2-x1)/mx;
    for j:=0 to my do
     begin
      initgraph(graphdriver,graphmode,'');
      ia:=0;ib:=round(sx*(xb-xa));ja:=maxy-round(-sy*ya);jb:=ja;
      line(ia,ja,ib,jb);
      ia:=round(-sx*xa);ib:=round(sx*(l-xa));
      ja:=jb;jb:=maxy-round(sy*(h-ya));line(ia,ja,ib,jb);
      ia:=ib;ib:=round(sx*(xb-xa));ja:=jb;line(ia,ja,ib,jb);
      for k:=0 to mx do
       begin
        x:=x1+k*(x2-x1)/mx;ia:=round(sx*(x-xa));ib:=ia;
        ja:=maxy-round(sy*(y1-ya));jb:=maxy-round(sy*(y2-ya));
        line(ia,ja,ib,jb);
       end;
      for k:=0 to my do
       begin
        y:=y1+k*(y2-y1)/my;ja:=maxy-round(sy*(y-ya));jb:=ja;
        ia:=round(sx*(x1-xa));ib:=round(sx*(x2-xa));
        line(ia,ja,ib,jb);
       end;
      if (h2<0.0) then
       begin
        ia:=0;ib:=round(sx*(xb-xa));ja:=maxy-round(sy*(h2-ya));jb:=ja;
        dots(ia,ja,ib,jb);
       end
      else
       begin
        ia:=0;ib:=round(sx*(h1*l/h-xa));ja:=maxy-round(sy*(h1-ya));jb:=ja;
        dots(ia,ja,ib,jb);
        ia:=round(sx*(h2*l/h-xa));ib:=round(sx*(xb-xa));
        ja:=maxy-round(sy*(h2-ya));jb:=ja;dots(ia,ja,ib,jb);
       end;
      yc:=y1+j*(y2-y1)/my;r:=yc-y0;
      p1:=1+h*h/(l*l);p2:=-2*xc-2*yc*h/l;p3:=xc*xc+yc*yc-r*r;
      xr:=(-p2+sqrt(p2*p2-4.0*p1*p3))/(2*p1);yr:=xr*h/l;
      if (xr>l) then
       begin
        yr:=h;xr:=xc+sqrt(r*r-(yc-yr)*(yc-yr));
       end;
      xl:=(-p2-sqrt(p2*p2-4.0*p1*p3))/(2*p1);yl:=xl*h/l;
      if (xl<0) then
       begin
        yl:=0.0;xl:=xc-sqrt(r*r-(yc-yl)*(yc-yl));
       end;
      if (yl<h1) then
       begin
        yl:=h1;xl:=xc-sqrt(r*r-(yc-yl)*(yc-yl));
       end;
      dx:=(xr-xl)/nx;
      ia:=round(sx*(xc-xa));ja:=maxy-round(sy*(yc-ya));
      x:=xl;yb:=yc-sqrt(r*r-(x-xc)*(x-xc));
      ib:=round(sx*(x-xa));jb:=maxy-round(sy*(yb-ya));
      line(ia,ja,ib,jb);
      for k:=1 to nx do
       begin
        x:=xl+(k-0.5)*dx;yb:=yc-sqrt(r*r-(x-xc)*(x-xc));
        ia:=ib;ja:=jb;ib:=round(sx*(x-xa));
        jb:=maxy-round(sy*(yb-ya));line(ia,ja,ib,jb);
        ta:=(x-xc)/(yc-yb);co:=sqrt(1.0/(1.0+ta*ta));si:=co*ta;
        yt:=x*h/l;if (x<0.0) then yt:=0.0;if (x>l) then yt:=h;
        if (h2<0.0) then yn:=h2 else
         begin
          x3:=h1*l/h;x4:=h2*l/h;
          yn:=x*h/l;
          if (x<x3) then yn:=h1;if (x>x4) then yn:=h2;
         end;
        p[k]:=0.0;if (yn>yb) then p[k]:=gw*(yn-yb);
        s[k]:=0.0;if (yb>0.0) then
         begin
          if (yn>yb) then
           begin
            if (yt>yn) then s[k]:=gn1*(yn-yb)+gd1*(yt-yn)
            else
             begin
              if (yt>yb) then s[k]:=gn1*(yt-yb)+gw*(yn-yt)
              else s[k]:=gw*(yn-yb);
             end;
           end
          else if (yt>yb) then s[k]:=gd1*(yt-yb);
         end
        else
         begin
          if (yn>0) then s[k]:=gn2*(0.0-yb) else
           begin
            if (yn>yb) then s[k]:=gn2*(yn-yb)+gd2*(0.0-yn)
            else s[k]:=gd2*(0.0-yb);
           end;
          if ((yt=0.0) and (yn>0)) then s[k]:=s[k]+gw*yn;
          if (yt>0.0) then
           begin
            if (yn>0.0) then
             begin
              if (yt>yn) then s[k]:=s[k]+gn1*yn+gd1*(yt-yn)
              else s[k]:=s[k]+gn1*yt+gw*(yn-yt);
             end
            else s[k]:=s[k]+gd1*yt;
           end;
         end;
       end;
      x:=xr;yb:=yc-sqrt(r*r-(x-xc)*(x-xc));
      ia:=ib;ja:=jb;ib:=round(sx*(x-xa));
      jb:=maxy-round(sy*(yb-ya));line(ia,ja,ib,jb);
      ia:=round(sx*(xc-xa));ja:=maxy-round(sy*(yc-ya));
      line(ia,ja,ib,jb);
      f:=1.0;e:=1.0;
      while (e>0.001) do
       begin
        a:=0.0;b:=0.0;
        for k:=1 to nx do
         begin
          x:=xl+(k-0.5)*dx;
          yb:=yc-sqrt(r*r-(x-xc)*(x-xc));
          ta:=(x-xc)/(yc-yb);co:=sqrt(1.0/(1.0+ta*ta));si:=co*ta;
          yt:=x*h/l;if (x<0.0) then yt:=0.0;if (x>l) then yt:=h;
          cc:=0.0;ph:=0.0;
          if (yb<0.0) then begin cc:=cc2;ph:=phi2;end
          else if (yb<yt) then begin cc:=cc1;ph:=phi1;end;
          tb:=sin(0.5*ph-0.25*pi)/cos(0.5*ph-0.25*pi);
          if (ta<tb) then ta:=tb;co:=sqrt(1.0/(1.0+ta*ta));
          ph:=pi*ph/180.0;tf:=sin(ph)/cos(ph);
          a:=a+(cc+(s[k]-p[k])*tf)/(co*(1+ta*tf/f));
          b:=b+s[k]*si;
         end;
        fa:=a/b;e:=abs(f-fa);f:=fa;
       end;
      ff[i,j]:=f;closegraph;
     end;
   end;
 end;
procedure output;
 var
  i,j:integer;
 begin
  title;
  writeln;writeln('     y    |    Stability factors :');
  writeln('          |');
  for j:=0 to my do
   begin
    write(y2-j*(y2-y1)/my:9:3,' |');
    for i:=0 to mx do write(ff[i,my-j]:9:3);writeln;
   end;
  writeln('    ----------------------------------------------------');
  write('     x    =');for i:=0 to mx do write(x1+i*(x2-x1)/mx:9:3);
  writeln;writeln;
 end;
begin
 graphinitialize;
 input;
 stability;
 output;
end.