Script MATLAB untuk Prediksi Elevasi Pasang Surut

Bongkar-bongkar harddisk dan laptop lama, ternyata saya menemukan sebuah script yang dibuat tahun 2011. Script ini dibuat untuk memprediksi elevasi pasang surut berdasarkan data amplitudo dan fasa dari 9 komponen pasang surut utama yang banyak digunakan di Indonesia. Daripada script ini hanya disimpan dan jamuran, berikut saya tuliskan saja di blog OSEANOGRAFI yang sudah jablay ini, siapa tahu bermanfaat buat mereka yang tidak sengaja berkunjung atau menemukan blog ini melalui mesin pencari.

Oh iya, di dalam script ini saya memanfaatkan code-code yang dibangun oleh Egbert dan Erofeeva (2002) dalam Oregon Tide Prediction Software (OTPS).

Enjoy!

% program prediksi elevasi pasang surut
% dengan input amplitudo (meter) dan fasa (derajat)
% untuk 9 komponen utama pasang surut
% m2,s2,k1,o1,n2,p1,k2
% agus setiawan - maret 2011

clear all

pi=3.14159265358979;
rad=pi/180.0;
ncmx=29;
nc=9; % jumlah komponen pasut untuk prediksi
ind=[1,2,3,4,5,6,7,21,22]; % nomor indeksnya dalam constid

% constid
% 'm2 ','s2 ','k1 ','o1 ',
% 'n2 ','p1 ','k2 ','q1 ',
% '2n2 ','mu2 ','nu2 ','l2 ',
% 't2 ','j1 ','m1 ','oo1 ',
% 'rho1','mf ','mm ','ssa ',
% 'm4 ','ms4 ','mn4 ','m6 ',
% 'm8 ','mk3 ','s6 ','2sm2',
% '2mk3'

% example data di stasiun tg. priok

% MSL
s0 = 0.0; % in meter

% amplitudo dan fasa
ampi(1) = 0.074; phasi(1) = 174.10; % m2
ampi(2) = 0.025; phasi(2) = 117.10; % s2
ampi(3) = 0.262; phasi(3) = 32.70; % k1
ampi(4) = 0.125; phasi(4) = 16.80; % o1
ampi(5) = 0.020; phasi(5) = 109.00; % n2
ampi(6) = 0.072; phasi(6) = 21.60; % p1
ampi(7) = 0.009; phasi(7) = 56.00; % k2
ampi(8) = 0.00; phasi(8) = 0.00; % m4
ampi(9) = 0.00; phasi(9) = 0.00; %ms4

% convert x and y components
for i=1:nc,
xx(i)=ampi(i)*cos(phasi(i)*rad);
yy(i)=ampi(i)*sin(phasi(i)*rad);
end

omega_d=[1.405189e-04,1.454441e-04,7.292117e-05,6.759774e-05,...
1.378797e-04,7.252295e-05,1.458423e-04,6.495854e-05,...
1.352405e-04,1.355937e-04,1.382329e-04,1.431581e-04,...
1.452450e-04,7.556036e-05,7.028195e-05,7.824458e-05,...
6.531174e-05,0.053234e-04,0.026392e-04,0.003982e-04,...
2.810377e-04,2.859630e-04,2.783984e-04,4.215566e-04,...
5.620755e-04,2.134402e-04,4.363323e-04,1.503693e-04,...
2.081166e-04];

phase_mkB=[1.731557546,0.000000000,0.173003674,1.558553872,...
6.050721243,6.110181633,3.487600001,5.877717569,...
4.086699633,3.463115091,5.427136701,0.553986502,...
0.052841931,2.137025284,2.436575100,1.929046130,...
5.254133027,1.756042456,1.964021610,3.487600001,...
3.463115091,1.731557546,1.499093481,5.194672637,...
6.926230184,1.904561220,0.000000000,4.551627762,...
3.809122439];

index=[30,35,19,12,27,17,37,10,25,26,28,33,34,...
23,14,24,11,5,3,2,45,46,44,50,0,42,51,40,0];

% awal tanggal peramalan
iyyy = 2010; % year
mm = 3; % month
id = 1; % day
hh = 0; % hour
mi = 0; % minute
ss = 0; % second
tgl0=datenum(iyyy,mm,id,hh,mi,ss);

% convert date to modified julian day
dpm(1) = 31; % jan
dpm(2) = 28; % feb
dpm(3) = 31; % mar
dpm(4) = 30; % apr
dpm(5) = 31; % may
dpm(6) = 30; % jun
dpm(7) = 31; % jul
dpm(8) = 31; % aug
dpm(9) = 30; % sep
dpm(10)= 31; % oct
dpm(11)= 30; % nov
dpm(12)= 31; % dec

mjd = 0;
% no earlier dates
if(iyyy<1858),iyyy=1858; end if(iyyy==1858 & mm>11),mm=11;end
if(iyyy==1858 & mm==11 & id>17),id=17;end

days=0;
for i=1:mm-1,
   days=days+dpm(i);
   if(i==2 & fix(iyyy/4)*4==iyyy),days=days+1;end
end
days=days+id-321;

% leap day correction
for k=1900:100:iyyy,
   if(iyyy==k & mm>2),days=days-1;end
end

for k=2000:400:iyyy,
   if(iyyy==k & mm>2),days=days+1;end
end

% each 4th year is leap year
nleap=fix((iyyy-1-1860)*0.25);
if(iyyy>1860),nleap=nleap+1;end

% except
for k=1900:100:iyyy-1,
   if(iyyy>k),nleap=nleap-1;end
   if(iyyy==k & mm>2),days=days-1;end
end

% but each in the row 2000:400:... is leap year again
for k=2000:400:iyyy-1,
   if(iyyy>k),nleap=nleap+1;end
   if(iyyy==k & mm>2),days=days+1;end
end
mjd=365*(iyyy-1858)+nleap+days;
time_mjd=mjd+hh/24.0+mi/(24.0*60.0)+ss/(24.0*60.0*60.0);

ntime=31*24*60; % lama prediksi nday*24hr*60 (in minute)
for k=1:1:ntime,
   tgl(k)=tgl0+k/1440;

   % astrol
   T=((k-1)/(24*60)+time_mjd)-51544.4993;
   circle=360.0;

   SHPN(1)=218.3164+13.17639648*T;
   SHPN(2)=280.4661+0.98564736*T;
   SHPN(3)=83.3535+0.11140353*T;
   SHPN(4)=125.0445-0.05295377*T;

   SHPN(1)=mod(SHPN(1),circle);
   SHPN(2)=mod(SHPN(2),circle);
   SHPN(3)=mod(SHPN(3),circle);
   SHPN(4)=mod(SHPN(4),circle);

   if(SHPN(1)<0.0),SHPN(1)=SHPN(1)+circle;end
   if(SHPN(2)<0.0),SHPN(2)=SHPN(2)+circle;end
   if(SHPN(3)<0.0),SHPN(3)=SHPN(3)+circle;end
   if(SHPN(4)<0.0),SHPN(4)=SHPN(4)+circle;end s=SHPN(1); 
   h=SHPN(2); 
   p=SHPN(3); 
   omega=SHPN(4); 
   pp=282.94; 

   % nodal 
   hour=(((k-1)/(24.0*60)+time_mjd)-fix(((k-1)/(24.0*60)+time_mjd)))*(24.0*60); 
   t1=15.0*hour; t2=30.0*hour; 
   arg(1)=h-pp; % Sa 
   arg(2)=2.0*h; % Ssa 
   arg(3)=s-p; % Mm 
   arg(4)=2.0*s-2.0*h; % MSf 
   arg(5)=2.0*s; % Mf 
   arg(6)=3.0*s-p; % Mt 
   arg(7)=t1-5.0*s+3.0*h+p-90.0; % alpha1 
   arg(8)=t1-4.0*s+h+2.0*p-90.0; % 2Q1 
   arg(9)=t1-4.0*s+3.0*h-90.0; % sigma1 
   arg(10)=t1-3.0*s+h+p-90.0; % q1 
   arg(11)=t1-3.0*s+3.0*h-p-90.0; % rho1 
   arg(12)=t1-2.0*s+h-90.0; % o1 
   arg(13)=t1-2.0*s+3.0*h+90.0; % tau1 
   arg(14)=t1-s+h+90.0; % M1 
   arg(15)=t1-s+3.0*h-p+90.0; % chi1 
   arg(16)=t1-2.0*h+pp-90.0; % pi1 
   arg(17)=t1-h-90.0; % p1 
   arg(18)=t1+90.0; % s1 
   arg(19)=t1+h+90.0; % k1 
   arg(20)=t1+2.0*h-pp+90.0; % psi1 
   arg(21)=t1+3.0*h+90.0; % phi1 
   arg(22)=t1+s-h+p+90.0; % theta1 
   arg(23)=t1+s+h-p+90.0; % J1 
   arg(24)=t1+2.0*s+h+90.0; % OO1 
   arg(25)=t2-4.0*s+2.0*h+2.0*p; % 2N2 
   arg(26)=t2-4.0*s+4.0*h; % mu2 
   arg(27)=t2-3.0*s+2.0*h+p; % n2 
   arg(28)=t2-3.0*s+4.0*h-p; % nu2 
   arg(29)=t2-2.0*s+h+pp; % M2a 
   arg(30)=t2-2.0*s+2.0*h; % M2 
   arg(31)=t2-2.0*s+3.0*h-pp; % M2b 
   arg(32)=t2-s+p+180.0; % lambda2 
   arg(33)=t2-s+2.0*h-p+180.0; % L2 
   arg(34)=t2-h+pp; % t2 
   arg(35)=t2; % 
   S2 arg(36)=t2+h-pp+180.0; % R2 
   arg(37)=t2+2.0*h; % K2 
   arg(38)=t2+s+2.0*h-pp; % eta2 
   arg(39)=t2-5.0*s+4.0*h+p; % MNS2 
   arg(40)=t2+2.0*s-2.0*h; % 2SM2 
   arg(41)=1.5*arg(30); % M3 
   arg(42)=arg(19)+arg(30); % MK3 
   arg(43)=3.0*t1; % S3 
   arg(44)=arg(27)+arg(30); % MN4 
   arg(45)=2.0*arg(30); % M4 
   arg(46)=arg(30)+arg(35); % MS4 
   arg(47)=arg(30)+arg(37); % MK4 
   arg(48)=4.0*t1; % S4 
   arg(49)=5.0*t1; % S5 
   arg(50)=3.0*arg(30); % M6 
   arg(51)=3.0*t2; % S6 
   arg(52)=7.0*t1; % S7 
   arg(53)=4.0*t2; 

   sinn=sin(omega*rad); 
   cosn=cos(omega*rad); 
   sin2n=sin(2.0*omega*rad); 
   cos2n=cos(2.0*omega*rad); 
   sin3n=sin(3.0*omega*rad); 

   f(1)=1.0; % Sa 
   f(2)=2.0; % Ssa 
   f(3)=1.0-0.130*cosn; % Mm 
   f(4)=1.0; % MSf 
   f(5)=1.043+0.414*cosn; % Mf 
   f(6)=sqrt((1.0+.203*cosn+.040*cos2n)^2+... 
        (.203*sinn+.040*sin2n)^2); % Mt 
   f(7)=1.0; % alpha1 
   f(8)=sqrt((1.+.188*cosn)^2+(.188*sinn)^2); % 2Q1 
   f(9)=f(8); % sigma1 
   f(10)=f(8); % q1 
   f(11)=f(8); % rho1 
   f(12)=sqrt((1.0+0.189*cosn-0.0058*cos2n)^2+...
         (0.189*sinn-0.0058*sin2n)^2); % O1 
   f(13)=1.0; % tau1 
   tmp1=1.36*cos(p*rad)+.267*cos((p-omega)*rad); % Ray's 
   tmp2= 0.64*sin(p*rad)+.135*sin((p-omega)*rad); 
   f(14)=sqrt(tmp1^2 + tmp2^2); % M1 
   f(15)=sqrt((1.+.221*cosn)^2+(.221*sinn)^2); % chi1 
   f(16)=1.0; % pi1 
   f(17)=1.0; % P1 
   f(18)=1.0; % S1 
   f(19)=sqrt((1.+.1158*cosn-.0029*cos2n)^2+...
         (.1554*sinn-.0029*sin2n)^2); % K1 
   f(20)=1.0; % psi1 
   f(21)=1.0; % phi1 
   f(22)=1.0; % theta1 
   f(23)=sqrt((1.+.169*cosn)^2+(.227*sinn)^2); % J1 
   f(24)=sqrt((1.0+0.640*cosn+0.134*cos2n)^2+...
         (0.640*sinn+0.134*sin2n)^2); % OO1 
   f(25)=sqrt((1.-.03731*cosn+.00052*cos2n)^2+...
         (.03731*sinn-.00052*sin2n)^2); % 2N2 
   f(26)=f(25); % mu2 
   f(27)=f(25); % N2 
   f(28)=f(25); % nu2 
   f(29)=1.0; % M2a 
   f(30)=f(25); % M2 
   f(31)=1.0; % M2b 
   f(32)=1.0; % lambda2 
   temp1=1.-0.25*cos(2.0*p*rad)...
        -0.11*cos((2.0*p-omega)*rad)-0.04*cosn; 
   temp2=0.25*sin(2.0*p)+0.11*sin((2.0*p-omega)*rad)...
        +0.04*sinn; 
   f(33)=sqrt(temp1^2+temp2^2); % L2 
   f(34)=1.0; % t2 
   f(35)=1.0; % S2 
   f(36)=1.0; % R2 
   f(37)=sqrt((1.+.2852*cosn+.0324*cos2n)^2+...
        (.3108*sinn+.0324*sin2n)^2); % K2 
   f(38)=sqrt((1.+.436*cosn)^2+(.436*sinn)^2); % eta2 
   f(39)=f(30)^2; % MNS2 
   f(40)=f(30); % 2SM2 
   f(41)=1.0; % M3 
   f(42)=f(19)*f(30); % MK3 
   f(43)=1.0; % S3 
   f(44)=f(30)^2; % MN4 
   f(45)=f(44); % M4 
   f(46)=f(44); % MS4 
   f(47)=f(30)*f(37); % MK4 
   f(48)=1.0; % S4 
   f(49)=1.0; % S5 
   f(50)=f(30)^3; % M6 
   f(51)=1.0; % S6 
   f(52)=1.0; % S7 
   f(53)=1.0; % S8 
   
   u(1)=0.0; % Sa 
   u(2)=0.0; % Ssa 
   u(3)=0.0; % Mm 
   u(4)=0.0; % MSf 
   u(5)=-23.7*sinn+2.7*sin2n-0.4*sin3n; % Mf 
   u(6)=atan(-(.203*sinn+.040*sin2n)/...
       (1.0+.203*cosn+.040*cos2n))/rad; % Mt 
   u(7)=0.0; % alpha1 
   u(8)=atan(.189*sinn/(1.+.189*cosn))/rad; % 2Q1 
   u(9)=u(8); % sigma1 
   u(10)=u(8); % q1 
   u(11)=u(8); % rho1 
   u(12)=10.8*sinn-1.3*sin2n+0.2*sin3n; % O1 
   u(13)=0.0; % tau1 
   u(14)=atan2(tmp2,tmp1)/rad; % M1 
   u(15)=atan(-.221*sinn/(1.+.221*cosn))/rad; % chi1 
   u(16)=0.0; % pi1 
   u(17)=0.0; % P1 
   u(18)=0.0; % S1 
   u(19)=atan((-.1554*sinn+.0029*sin2n)/...
        (1.+.1158*cosn-.0029*cos2n))/rad; % K1 
   u(20)=0.0; % psi1 
   u(21)=0.0; % phi1 
   u(22)=0.0; % theta1 
   u(23)=atan(-.227*sinn/(1.+.169*cosn))/rad; % J1 
   u(24)=atan(-(.640*sinn+.134*sin2n)/...
        (1.+.640*cosn+.134*cos2n))/rad; % OO1 
   u(25)=atan((-.03731*sinn+.00052*sin2n)/...
        (1.-.03731*cosn+.00052*cos2n))/rad; % 2N2 
   u(26)=u(25); % mu2 
   u(27)=u(25); % N2 
   u(28)=u(25); % nu2 
   u(29)=0.0; % M2a 
   u(30)=u(25); % M2 
   u(31)=0.0; % M2b 
   u(32)=0.0; % lambda2 
   u(33)=atan(-temp2/temp1)/rad; % L2 
   u(34)=0.0; % t2 
   u(35)=0.0; % S2 
   u(36)=0.0; % R2 
   u(37)=atan(-(.3108*sinn+.0324*sin2n)/...
        (1.+.2852*cosn+.0324*cos2n))/rad; % K2 
   u(38)=atan(-.436*sinn/(1.+.436*cosn))/rad; % eta2 
   u(39)=u(30)*2.0; % MNS2 
   u(40)=u(30); % 2SM2 
   u(41)=1.50*u(30); % M3 
   u(42)=u(30)+u(19); % MK3 
   u(43)=0.0; % S3 
   u(44)=u(30)*2.0; % MN4 
   u(45)=u(44); % M4 
   u(46)=u(30); % MS4 
   u(47)=u(30)+u(37); % MK4 
   u(48)=0.0; % S4 
   u(49)=0.0; % S5 
   u(50)=u(30)*3.0; % M6 
   u(51)=0.0; % S6 
   u(52)=0.0; % S7 
   u(53)=0.0; % S8 

   pu=zeros(1,ncmx); 
   pf=ones(1,ncmx); 

   for i=1:ncmx, 
      if(index(i)>0),
         pu(i)=u(index(i))*rad;
         pf(i)=f(index(i));
      end
   end

   % end of nodal

   time=(((k-1)/(24.0*60)+time_mjd)-48622.0)*86400;

   % make_a

   for j=1:ncmx,
      omegi(j)=omega_d(j);
      phase(j)=phase_mkB(j);
   end

   for j=1:nc,
      i=ind(j);
      if(i~=0),
         xc(j)=pf(i)*cos(omegi(i)*time+phase(i)+pu(i));
         yc(j)=pf(i)*sin(omegi(i)*time+phase(i)+pu(i));
      end
   end

   sum=0.0;
   for i=1:nc,
      sum=sum+xx(i)*xc(i)+yy(i)*yc(i);
   end
   zpred(k)=sum;
end

% gambarkan hasil prediksi
dateFormat='dd/mm/yyyy';
figure('Color','w');
plot(tgl,zpred);
grid on
set(gca, 'YLim', [-0.5, 0.5]);
set(gca, 'XLim', [tgl0 tgl0+15]);
title('Tidal Elevation Prediction');
ylabel('Elevation (m)');
datetick('x',dateFormat)

 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s