******************************************************
*          ストックホルムの潮汐の解析 -Fortran-      *
*                        - Goshi Fujimoto -          *
******************************************************

      program stock

      parameter(yrn=2003-1889+1)
      parameter(n=yrn*12)
      parameter( ID0=18890101, ND=1380, RND=ND )
      
      real tyo(yrn,12)
      real date(n), tyo2(n), reg(n), date2(n)
      integer iws, yr, i
      integer year(yrn)
      real datetotal, tyo2total, datetyo2, date2jou, a, k


*Read of data

      open(10,file='./tyoui.dat',form='formatted')

      do yr=1,yrn
         read(10,*) year(yr), (tyo(yr,mon), mon=1,12) !データ読み込み(年回し)。1行目の年とその後の年月の潮位。行列みたいな形で読み込んでる。
      end do

!読み込んだデータを２次元にする。そのために新しい変数iを導入。

      do yr=1, yrn
       do mon=1, 12
          i = (yr-1)*12 + mon
          tyo2(i) = tyo(yr,mon) !2次元にした潮位
          date(i) = real(i) !2次元にしたデータ数or日
       end do
      end do


*Regression line by least square
      
      datetotal = 0.
      tyo2total = 0.
      datetyo2 = 0.
      date2jou = 0.
      k = 0. !足したり引いたりするときは初期値を設定した方がいい。
      do i=1, n
         if (tyo2(i).ne.-999.) then !欠損値999でないときのみ計算
            datetotal = datetotal + date(i) !日数総数。割ることで平均日数
            tyo2total = tyo2total + tyo2(i) !潮位値総数。割る平均潮位
            datetyo2 = datetyo2 + date(i)*tyo2(i) !日と潮位値の掛算の総足し
            date2jou = date2jou + date(i)*date(i) !日^2乗の総数
            k = k + 1
         endif
      enddo

      a = ((datetyo2/k) - ((datetotal/k) * (tyo2total/k)))/
     &        ((date2jou/k) - (datetotal/k)*(datetotal/k)) !回帰線の傾き

      do i=1, n
         reg(i) = a*date(i) + ((tyo2total/k) - (a*datetotal/k)) !回帰直線
         date2(i) = (date(i)/12) + 1889 !x軸を年単位にするために1380数をもつdateを115年という数にここで変換
      enddo

c      write(10, rec=1) reg

*Configuration of DCL
      write(*,*) 'Workstation ID?'
      call SGPWSN                    
      read (*,*) iws                 !デバイス番号を入力
      call SGOPN( iws )                !デバイスをオープンする 
      
      call UZFACT( 0.8 ) !字やメモリの倍率
      call UZLSET( 'LOFFSET', .TRUE. ) !軸を二つ付ける準備。他に使い道あり。
      call SLMGN( 0.05, 0.05, 0.05, 0.05 ) !マージンの領域
c      call SLSTTL( 'Tide level at Stockholm', 'T', 0., -1., 0.03, 1 )
      call SLSTTL( 'Goshi Fujimoto', 'B', 0., 0., 0.02, 2 )
      call SLSTTL( '#date #time', 'B', 1., -1., 0.01, 4 )
      call GLLSET( 'lmiss', .true. )    !欠損値処理をONにする
      call GLRSET( 'rmiss', -999. )    !欠損値を -999.0に設定する

      call GRFRM !新しい作画領域を設定
      call GRSWND( 1889., 2004., 6400., 8000. ) !ウインドウパラメータの設定
      call GRSVPT( 0.2, 0.85, 0.2, 0.8 ) !ビューポートの設定
      call GRSTRN( 1 ) !正規化変換の関数番号
      call GRSTRF !正規化変換の確定

c      call UCXACL( 'B', ID0, ND ) !x軸の設定(時間軸)
c      call UCXACL( 'T', ID0, ND )
c      call UXSAXS( 'B' ) !同じ軸に何本も軸を
      call UXAXDV( 'T', 5., 20. ) !x軸の設定(date2に対する軸)
      call UXAXDV( 'B', 5., 20. )

      call UYAXDV( 'L', 50., 200. ) !y軸の設定
      call UYAXDV( 'R', 50., 200. )

      call UXSTTL( 'B', '[Year]', 1. ) !軸の単位
      call UXMTTL( 'B', 'Date', 0. ) !軸のラベル
      call UYSTTL( 'L', '[m]', 1. )
      call UYMTTL( 'L', 'Tide Level', 0. )

*     ..... make figure .....
c      call UUSLNT( 1 ) !1:実線 2:破線 3:点線
c      call UUSLNI( 23 ) !dclclr参照
      call UULINZ( n, date2, tyo2, 1, 111 ) !date2で作図
c      call UUSLNT( 1 ) !1:実線 2:破線 3:点線
c      call UUSLNI( 403 ) !dclclr参照
      call UULINZ( n, date2, reg, 1, 202 )

*     ..... legend .....
*     Max
      call SGSLNI( 1 )
      call SGSLNI( 113 )
      call SGLNU( 1960., 7800., 1970., 7800. )
*     Min
      call SGSLNI( 1 )
      call SGSLNI( 202 )
      call SGLNU( 1960., 7700., 1970., 7700. )
*     font
      call SGSTXS( 0.015 )
      call SGSTXC( -1 )
      call SGTXU( 1973., 7800., 'Tide Level' )
      call SGTXU( 1973., 7700., 'Reg. Line' )

      call GRCLS !figure finish

      stop
      end 
