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)