matlab实现海洋要素计算代码
基于MATLAB的潮流计算源程序代码(优.选)
%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算**********clearclcload E:\data\IEEE014_Node.txtNode=IEEE014_Node;weishu=size(Node);nnum=weishu(1,1); %节点总数load E:\data\IEEE014_Branch.txtbranch=IEEE014_Branch;bwei=size(branch);bnum=bwei(1,1); %支路总数Y=(zeros(nnum));Sj=100;%********************************节点导纳矩阵*******************************for m=1:bnum;s=branch(m,1); %首节点e=branch(m,2); %末节点R=branch(m,3); %支路电阻X=branch(m,4); %支路电抗B=branch(m,5); %支路对地电纳k=branch(m,6);if k==0 %无变压器支路情形Y(s,e)=-1/(R+j*X); %互导纳Y(e,s)=Y(s,e);endif k~=0 %有变压器支路情形Y(s,e)=-(1/((R+j*X)*k));Y(e,s)=Y(s,e);Y(s,s)=-(1-k)/((R+j*X)*k^2);Y(e,e)=-(k-1)/((R+j*X)*k); %对地导纳endY(s,s)=Y(s,s)-j*B/2;Y(e,e)=Y(e,e)-j*B/2; %自导纳的计算情形endfor t=1:nnum;Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13);%求支路自导纳endG=real(Y); %电导B=imag(Y); %电纳%******************节点分类************************************* *pq=0; pv=0; blancenode=0;pqnode=zeros(1,nnum);pvnode=zeros(1,nnum);for m=1:nnum;if Node(m,2)==3blancenode=m; %平衡节点编号else if Node(m,2)==0pq=pq+1;pqnode(1,pq)=m; %PQ 节点编号else if Node(m,2)==2pv=pv+1;pvnode(1,pv)=m; %PV 节点编号endendendend%*****************************设置电压初值********************************** Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化for n=1:nnumUoriginal(1,n)=Node(n,9); %对各点电压赋初值if Node(n,9)==0;Uoriginal(1,n)=1; %该节点为非PV节点时，将电压值赋为1endendPresion=input('请输入误差精度要求：Presion=');disp('该电力系统节点数:');disp(nnum);xiumax=0.1;counter=0;while xiumax>Presion%****************************计算不平衡量***********************************e=real(Uoriginal); %取初始电压的实部f=imag(Uoriginal); %取初始电压的虚部deta=zeros(2*pq+2*pv,1); %构造储存功率变化量的列矩阵n=1;for m=1:pq;Pi=0;Qi=0;for t=1:nnum;Pi=Pi+e(1,pqnode(1,m))*(G(pqnode(1,m),t)*e (1,t)-B(pqnode(1,m),t)*f(1,t))+f(1,pqnode(1,m ))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t) *e(1,t));%计算该PQ节点的负荷有功Qi=Qi+f(1,pqnode(1,m))*(G(pqnode(1,m),t)*e (1,t)-B(pqnode(1,m),t)*f(1,t))-e(1,pqnode(1,m ))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t) *e(1,t));%计算该PQ节点的负荷无功endS1(1,pqnode(1,m))=Pi+j*Qi;P=(Node(pqnode(1,m),7)-Node(pqnode(1,m), 5))/Sj-Pi;%计算该PQ节点的实际有功功率deta(n,1)=P;%在该列向量中储存有功功率n=n+1;Q=(Node(pqnode(1,m),8)-Node(pqnode(1,m),6))/Sj-Qi;%计算该PQ节点的实际无功功率deta(n,1)=Q;%在该列向量中储存无功功率n=n+1;endfor m=1:pv;Pv=0; Qv=0;for t=1:nnum;Pv=Pv+e(1,pvnode(1,m))*(G(pvnode(1,m),t)* e(1,t)-B(pvnode(1,m),t)*f(1,t))+f(1,pvnode(1, m))*(G(pvnode(1,m),t)*f(1,t)+B(pvnode(1,m), t)*e(1,t));%计算该PV节点的负荷有功Ui=e(1,pvnode(1,m))^2+f(1,pvnode(1,m))^2; %计算该节点的负荷电压值Qv=Qv+f(1,pqnode(1,m))*(G(pqnode(1,m),t)* e(1,t)-B(pqnode(1,m),t)*f(1,t))-e(1,pqnode(1, m))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m), t)*e(1,t));endS1(1,pvnode(1,m))=Pv+j*Qv;P=(Node(pvnode(1,m),7)-Node(pvnode(1,m), 5))/Sj-Pv; %计算该节点的实际有功功率deta(n,1)=P; %储存该有功功率n=n+1;U=Node(pvnode(1,m),3)^2-Ui; %计算电压变化量deta(n,1)=U; %储存该电压变化量n=n+1;enddeta;cerate=zeros(pq+pv,1);for k=1:pqcerate(k,1)=pqnode(1,k);endfor v=1:pvcerate(pq+v,1)=pvnode(1,v);end%******************************雅克比矩阵****************************** Jacob=ones(2*nnum-2);L=0;J=0;H=0;N=0; R=0;S=0;n=1;k=1;form=1:pq; %m表示雅克比矩阵中pq节点的行数for u=1:pq+pv; %u 表示雅克比矩阵中pq节点的列数t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的个数if pqnode(1,m)~=t %非对角元素的情况H=G(pqnode(1,m),t)*f(1,pqnode(1,m))-B(pqn ode(1,m),t)*e(1,pqnode(1,m));N=G(pqnode(1,m),t)*e(1,pqnode(1,m))+B(pq node(1,m),t)*f(1,pqnode(1,m));L=H;J=-N;elseifpqnode(1,m)==t %对角线元素时的情况I=0;for g=1:nnumI=Y(t,g)*Uoriginal(1,g)+I; %计算节点的注入电流endaii=real(I);bii=imag(I);H=-B(t,t)*e(1,pqnode(1,m))+G(t,t)*f(1,pqnode (1,m))+bii;N=G(t,t)*e(1,pqnode(1,m))+B(t,t)*f(1,pqnode (1,m))+aii;L=-B(t,t)*e(1,pqnode(1,m))+G(t,t)*f(1,pqnode (1,m))-bii; J=-G(t,t)*e(1,pqnode(1,m))-B(t,t)*f(1,pqnode( 1,m))+aii;endendJacob(n,k)=H;k=k+1;Jacob(n,k)=N;k=k-1;n=n+1;Jacob(n,k)=J;k=k+1;Jacob(n,k)=L;n=n-1; k=k+1; %按照雅克比矩阵的排列规则排列pq节点的雅克比元素endk=1; n=2*m+1; %将光标定位于下一个待排列PQ节点元素的第一个位置endn=2*pq+1; k=1; %定位于PV节点的第一个位置处for m=1:pv;for u=1:pq+pv;t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的位置if pvnode(1,m)~=t %非对角线元素情况H=G(pvnode(1,m),t)*f(1,pvnode(1,m))-B(pvno de(1,m),t)*e(1,pvnode(1,m));N=G(pvnode(1,m),t)*e(1,pvnode(1,m))+B(pvn ode(1,m),t)*f(1,pvnode(1,m));R=0; S=0;endifpvnode(1,m)==t %对角线元素情况I=0;for g=1:nnumI=Y(t,g)*Uoriginal(1,g)+I; %计算PV节点的注入电流endaii=real(I);bii=imag(I);H=-B(t,t)*e(1,pvnode(1,m))+G(t,t)*f(1,pvnode (1,m))+bii;N=G(t,t)*e(1,pvnode(1,m))+B(t,t)*f(1,pvnode( 1,m))+aii;R=2*f(1,pvnode(1,m));S=2*e(1,pvnode(1,m));endJacob(n,k)=H; k=k+1;Jacob(n,k)=N; k=k-1;n=n+1;Jacob(n,k)=R; k=k+1;Jacob(n,k)=S;n=n-1;k=k+1; %按照雅克比矩阵的排列规则排列PV节点的雅克比元素endk=1;n=n+2; %定位于下一个待排列PV节点的雅克比元素第一个位置end%*************************电压变化量化的计算与存储************************************ Detau=inv(Jacob)*deta; %构建电压的变化量的列向量f=zeros(1,nnum); %给电压实部赋初值0e=zeros(1,nnum); %给电压虚部赋初值0for p=1:(pq+pv);f(1,cerate(p,1))=j*Detau(2*p-1,1);%将电压变量的奇数行赋值给fe(1,cerate(p,1))=Detau(2*p,1); %将电压变量的偶数行赋值给eendt=e+f;xiumax=abs(Detau(1,1)); %将电压变化量的第一个元素赋值给最大允许误差for n=2:2*nnum-2;if abs(Detau(n,1))>xiumaxxiumax=abs(Detau(n,1)); %找出最大的电压误差endendUoriginal=Uoriginal+t; %迭代修正后的电压值counter=1+counter; %统计迭代次数enddisp('迭代次数counter:');disp(counter);%**************************平衡节点功率及显示**********************************m=blancenode;t=0;for n=1:nnum;t=t+(G(m,n)-j*B(m,n))*(real(Uoriginal(1,n))-j*i mag(Uoriginal(1,n)));endS1(1,m)=Uoriginal(1,m)*t;%**************************直角坐标下各节点电压及显示****************************U=zeros(1,nnum);for n=1:nnumUi(n,1)=Node(n,1);U(1,n)=real(Uoriginal(1,n))+i*imag(Uoriginal(1 ,n)); %将电压值由极坐标转化为直角坐标形式Ui(n,2)=U(1,n);Ui(n,3)=S1(1,n);enddisp('各节点电压直角坐标形式及节点注入功率:');disp(' 节点号节点电压值节点注入功率');disp(Ui);disp('修正电压的最大误差:')disp(xiumax);%**************************功率损耗************************************* **for m=1:bnum %支路功率及损耗startnode=branch(m,1);endnode=branch(m,2); %终止节点y=sum(Y,2);Sij=Uoriginal(1,startnode)*((conj(Uoriginal(1,s tartnode))*conj(y(startnode,1)))+(conj(Uorigi nal(1,startnode))-conj(Uoriginal(1,endnode))) *conj(-Y(startnode,endnode)));Sji=Uoriginal(1,endnode)*((conj(Uoriginal(1,e ndnode))*conj(y(endnode,1)))+(conj(Uoriginal (1,endnode))-conj(Uoriginal(1,startnode)))*co nj(-Y(endnode,startnode)));S(m,1)=startnode; %始节点S(m,2)=endnode; %末节点S(m,3)=Sij; %始节点流入的功率S(m,4)=Sji; %末节点流入的功率S(m,5)=Sij+Sji; %线路损耗enddisp('支路功率及损耗：');disp('节点号(I) 节点号(J) 支路功率（I-J）支路功率（J-I）线路损耗（delta_S)')disp(S); 最新文件---------------- 仅供参考--------------------已改成word文本--------------------- 方便更改。
基于VB和Surfer的海洋要素制图的批处理可视化系统
22科技资讯 SCIENCE & TECHNOLOGY INFORMATION信 息 技 术①基金项目：国家海洋局南海分局海洋科学技术局长基金资助(项目编号：1519)。
作者简介：钟煜宏(1985—)，男，汉族，广东河源人，本科，助理工程师，研究方向：海洋信息技术。
通讯作者：吴梅桂(1986—)，女，汉族，福建泉州人，硕士，助理工程师，研究方向：海洋放射化学、海洋信息技术，E- mail:rosemg1027@。
DOI：10.16661/ki.1672-3791.2018.23.022基于VB和Surfer的海洋要素制图的批处理可视化系统①钟煜宏 吴梅桂*(国家海洋局南海环境监测中心 广东广州 510300)摘 要：本文基于VB和Automation对象技术的使用程序控制Surfer自动绘图的方法，利用Surfer自带的Scripter脚本语言，可以实现Surfer面向对象编程语言的二次开发，实现海洋要素等值线图、分类图等不同类型图件的批量绘制,使用Visual Studio 2013开发工具，实现批处理可视化系统功能和界面，系统提供海洋数据处理功能，图件批处理等操作界面。
本系统可以显著提高海洋绘图工作者的工作效率，避免重复工作导致的误差。
关键词：Surfer Automation VB 等值线 批量绘图 海洋要素中图分类号：TP31 文献标识码：A 文章编号：1672-3791(2018)08(b)-0022-03海洋制图数据量巨大，如果手动一个个处理数据，画等值线，效率低下，甚至容易出错，导致模拟结果失真。
根据需要，本文比较几种处理方法，形成批量画等值线图的可视化系统。
Surfer是美国Golden Software公司研制开发的绘图软件,其不仅提供了丰富的网格化和插值方法,还具有强大的绘制等值线等矢量图能力，而且Sur fer的还提供了Automation技术,发现可以使用Visual Studio 2010的语言开发批量制图的客户端系统，此程序后台调用Surfer的Automation对象，轻松实现Surfer的强大图形绘制功能[1]，本文选择Surfer软件的Automation技术去实现水质模拟可视化系统，去处理批量的数据。
Matlab在水资源规划中的应用
Matlab在水资源系统规划中的应用缴锡云§1 Matlab初步§1.1 数值计算功能1 基本数学运算◆符号“%”后为注解；◆表达式后加分号—“；”，回车后不显示结果；◆表达式续行，可以用“…”连接；◆运算的优先级规则符合一般规则；◆Matlab中所使用的逗号、分号等必须为半角（在英文状态下输入）。
表1.1 代数运算表运算符号举例加法：x＋y + 1+2减法：x－y - 100-50乘法：x×y * 9.8*2除法：x÷y /或\ 65/5=5\65乘幂：x y^ 3^2%示例1—1：1+2+3ans=6%示例1—2：x=1+2+3x=6%示例1—3：x=2;y=3;z=x*y+2z=82 变量◆变量命名的规则（1）字母打头，可含数字、下划线；（2）区分字母大小写；（3）不超过19个字符。
◆变量的删除clear 变量名◆特殊变量表1.2 特殊变量表特殊变量变量说明ans 结果的缺省变量名pi 圆周率eps 系统的最小数inf（或Inf）无穷大i、j 虚数单位注：特殊变量可重新定义；表中，除了Inf以外，其余均小写有效%示例1—4：pians=3.1416pi=2pi=2clear pipi???Undefined function or variable‘pi’3 数学函数表1.3 常用函数表函数说明abs(x) 绝对值acos(x) 反余弦asin(x) 反正弦atan(x) 反正切cos(x) 余弦sin(x) 正弦tan(x) 正切exp(x) 指数函数log10(x) 常用对数表1.3 常用函数表（续）函数说明log(x) 自然对数sqrt(x) 开平方rem(x,y) 求x/y的余数fix(x) 接近0取整floor(x) 接近负无穷取整round(x) 四舍五入取整sign(x) 符号函数：x<0,sign(x)=-1;x=0,sign(x)=0;x>0,sign(x)=14 关系与逻辑运算在执行关系及逻辑运算时，MA TLAB 将输入的不为零的数值都视为真（True），而为零的数值则视为否（False）。
海洋要素计算(潮汐)
海洋要素计算作业之二——潮汐（威海2013年五月份）一．本次潮汐调和分析共选取了十三个分潮：MSf，Q1，O1，K1，P1，K2，N2，M2，S2，MK3，M4，MS4，M6为使您查看方便，将本次大作业的放在本文件夹各文件内，具体参考如下：1.原数据为：qd.dat;2.Fortran编程见该文件夹内：tide.f90文件；3.求各分潮调和常数H、g的值及其中间过程得到的各值见：qd_tide.dat文件；二．对比回报值和实测值：1. 回报1968年一月份的水位值见：huibao.dat；2. 用matlab绘制的潮汐过程曲线见：潮汐过程曲线.bmp3. 用给定的六个分潮求得的高潮和低潮发生的时刻及潮位值见—：gaodichao.dat;运行tide.f90后求得威海地区2013年5月份的平均潮差。
由图可知：由于只计算了一个月的潮汐数据，所以回报值和实测值相符的不是很好，如果计算一年的数据，应该会取得比较良好的结果。
三．程序%% 潮汐过程曲线图clear,clc%%huibao=load('G:\chaoxi\huibao.dat');% huibao=fread(fhuibao);shice=load('G:\chaoxi\qd.dat');% shice=fread(fshice);%huibao_y=zeros(1,12*62);%shice_y=zeros(1,12*62);huibao=double(huibao');huibao_y=double(huibao(:));%shice_y=reshape(shice',1,[])%for i=1:12;% for j=1:62% huibao_y(i)=huibao(i,j)% shice_y(i)=shice(i,j)%end%endshice=double(shice');shice_y=double(shice(:));x=linspace(1,31,length(huibao_y));plot(x,huibao_y,'r-')hold onplot(x,shice_y,'b-')title('威海(37°31′N ,122°08′E)2013年五月潮汐调和分析图') legend('回报值','实测值')xlabel('时间（2013年五月份）')ylabel('水位（m）')。
Matlab.NET混合编程在潮汐调和分析与预报中的应用
关键词 ： Ｍ ａ ｔ ｌ ａ ｂ ． Ｎ Ｅ Ｔ混编 ； Ｔ ＿ Ｔ Ｉ Ｄ Ｅ程序 包 ； Ｄ Ｌ Ｌ （ 动 态链 接 库 ） ； 潮 汐调 和 分 析 ； 潮位 预 报
合编程技术 ， 创 建潮汐调 和常数计算 Ｔ Ｈ Ｃ Ｃ ａ ｌ ｃ ｕ ｌ ａ ｔ ｅ ． ｄ ｌ ｌ 和潮 汐预报 Ｔ ｉ ｄ ｅ Ｆ ｏ ｒ ｅ ｃ ａ ｓ ｔ ． ｄ ｌ ｌ ， 在主程序 中调用对应 的函数 ， 实现 了操作 简单且具有 良好 用户界 面的潮 汐调 和分析与潮位预报 系统 ， 具有潮位数 据读取、 调和分析 和预 报、 数
（ １ ． Ｃ ｏ ｌ ｌ ｅ ｇ ｅ ｏ ｆ Ｇｅ ｏ ｍａ ｔ ｉ ｃ ｓ ，Ｓ ｈ ａ ｎ ｄ ｏ ｎ ｇ Ｕｎ ｉ ｖ ｅ ｒ ｓ ｉ ｔ ｙ ｏ ｆ Ｓ ｃ ｉ ｅ ｎ ｃ ｅ ａ ｎ ｄ Ｔ ｅ ｃ ｈ ｎ ｏ ｌ ｏ ｇ ｙ ， Ｑｉ ｎ ｇ ｄ ａ ｏ ２ ６ ６ ５ ９ ０ ，Ｃ ｈ ｉ ｎ ａ ；
ｉ ｎ ｔ ｅ ｒ ｆ ａ ｃ ｅ ａ ｎ ｄ ｔ ｉ ｄ ｅ ｄ ａ ｔ ａ ｒ ｅ ａ ｄ ｉ ｎ ｇ ，ｈ ａ ｒ ｍｏ ｎ ｉ ｃ ａ ｎ ｌｙ ａ ｓ ｉ ｓ ａ ｎ ｄ ｆ ｏ ｒ ｅ ｃ ａ ｓ ｔ ｉ ｎ ｇ ，ｄ ａ ｔ ａ ｒ ｅ ｎ ｄ ｅ ｉｎ ｒ ｇ， ｅ ｔ ｃ ．ｂ ｙ ｕ ｓ ｉ ｎ ｇ ｔ ｈ ｅ Ｍａ ｌ ｆ ａ ｂ ｏ ｒａ ｇ ｍ－
中图分类号 ： Ｐ ２ ０ ９
文献标 识码 ： Ａ
文章编 号 ： １ ６ ７ ２ — ５ ８ ６ ７ （ ２ ０ １ ７ ） ０ ５— ０ ０ ５ ９— ０ ３
人工鱼群算法matlab实现
function lhl_AFclc;clear all; close all;format longVisual = 2.5; %人工鱼的感知距离Step = 0.3; %人工鱼的移动最大步长N = 10; %人工鱼的数量Try_number = 50;%迭代的最大次数delta=0.618; %拥挤度因子a1 = -10; b1 = 10; a2 = -10; b2 = 10;d = [];%存储50个状态下的目标函数值;k = 0;m = 50;%迭代次数X1 = rand(N,1)*(b1-a1)+a1; %在-10~10之间，随机生成50个数；X2 = rand(N,1)*(b2-a2)+a2;X = [X1 X2];%X = ones(N,2);%for i = 1:N% X(i,1)=-10;% X(i,2)=10;%end% 人工鱼数量,两个状态变量X1和X2；%计算50个初始状态下的；for i = 1:Nwww = [X(i,1),X(i,2)];d(i) = maxf(www);end%公告牌用于记录人工鱼个体的历史最好状态[w,i] = max(d); % 求出初始状态下的最大值w和最大值的位置i；maxX = [X(i,1),X(i,2)]; % 初始公告板记录，最大值位置；maxY = w; % 初始化公告板记录，最大值；figurex = []; figurey = []; figurez = [];figurex(numel(figurex)+1) = maxX(1); % 将maxX（1）放入figurex中，figurey(numel(figurey)+1) = maxX(2); % numel返回数组或者向量中所含元素的总数，matlab数组下标默认是从1开始的figurez(numel(figurez)+1) = maxY;while(k<m)for i = 1:NXX = [X(i,1),X(i,2)]; %拿出其中一条鱼来看他的四种行为判断%%%%%%第一种行为：聚群行为：伙伴多且不挤，就向伙伴中心位置移动%群聚行为是伙伴的中心点，凸规划下，中心点一定还在约束内%群聚行为不是一种maxf（Xc）的比较，就是看伙伴位置nf1=0;Xc=0;label_swarm =0; %群聚行为发生标志for j = 1:NXX_1 = [X(j,1), X(j,2)];if (norm(XX_1-XX)<Visual) % norm函数求向量XXX-XX的范数，由于二维向量，2或者省略都可以nf1 = nf1+1;Xc = Xc+XX_1;endendXc=Xc-XX; %需要去除XX本身；nf1=nf1-1;Xc = Xc/nf1; %此时Xc表示XX感知范围其他伙伴的中心位置；if((maxf(Xc)/nf1 > delta*maxf(XX)) && (norm(Xc-XX)~=0))XXR1=rand*Step*(Xc-XX)/norm(Xc-XX);XXnext1=XX+XXR1;if(XXnext1(1) > b1)XXnext1(1) = b1;endif(XXnext1(1) < a1)XXnext1(1) = a1;endif(XXnext1(2) > b2)XXnext1(2) = b2;endif(XXnext1(2) < a2)XXnext1(2) = a2;endlabel_swarm =1;temp_y_XXnext1=maxf(XXnext1);elselabel_swarm =0;temp_y_XXnext1=-inf;end%%%%%%%%%%%%第二种行为：追尾行为：周围伙伴有最大值且附近不挤，向其伙伴方向移动%追尾行为追寻伙伴行为，还是在约束内temp_maxY = -inf; %按照理论来说这块应该初始化为-无穷小，label_follow =0;%追尾行为发生标记for j = 1:NXX_2 = [X(j,1),X(j,2)];if((norm(XX_2-XX)<Visual) && (maxf(XX_2)>temp_maxY))temp_maxX = XX_2;temp_maxY = maxf(XX_2);endendnf2=0;for j = 1:NXX_2 = [X(j,1),X(j,2)];if(norm(XX_2-temp_maxX)<Visual)nf2=nf2+1;endendnf2=nf2-1;%去掉他本身if((temp_maxY/nf2)>delta*maxf(XX) && (norm(temp_maxX-XX)~=0)) %附近有Yj最大的伙伴，并且不太拥挤XXR2=rand*Step*(temp_maxX-XX)/norm(temp_maxX-XX);%rand不是随机反向，是随机步长XXnext2 = XX+XXR2;if(XXnext2(1) > b1)XXnext2(1) = b1;endif(XXnext2(1) < a1)XXnext2(1) = a1;endif(XXnext2(2) > b2)XXnext2(2) = b2;endif(XXnext2(2) < a2)XXnext2(2) = a2;endlabel_follow =1;temp_y_XXnext2=maxf(XXnext2);elselabel_follow =0;temp_y_XXnext2=-inf;end%%%%%%%%%%%%第三种行为：觅食行为：与前两个行为不同，觅食和随机行为都是找附近的状态，而不是找附近的同伴%觅食和随机行为可能出现超出约束，所以，XX_3和XX_4是不一样的%觅食行为和群聚行为、追尾行为是不一样的，觅食行为是一种根据状态来判断的行为，群聚和追尾是根据伙伴来判断的行为label_prey =0; %判断觅食行为是否找到优于当前的状态for j = 1:Try_numberR1V=Visual*(-1+2*rand(2,1)');XX_3 = XX+R1V;if(XX_3(1) > b1) % 下面这四个是一套，如果超出约束条件，就选值为边界条件XX_3(1) = b1;endif(XX_3(1) < a1)XX_3(1) = a1;endif(XX_3(2) > b2)XX_3(2) = b2;endif(XX_3(2) < a2)XX_3(2) = a2;endif(maxf(XX)<maxf(XX_3))XXR3=rand*Step*(XX_3-XX)/norm(XX_3-XX);XXnext3 = XX+XXR3;if(XXnext3(1) > b1) % 下面这四个是一套，如果超出约束条件，就选值为边界条件XXnext3(1) = b1;endif(XXnext3(1) < a1)XXnext3(1) = a1;endif(XXnext3(2) > b2)XXnext3(2) = b2;endif(XXnext3(2) < a2)XXnext3(2) = a2;endlabel_prey =1;break;endendtemp_y_XXnext3=max(XXnext3);if(label_prey ==0)temp_y_XXnext3=-inf;end%%%%%%%%%%%%行为选择if((label_swarm==0) && (label_follow==0) && (label_prey ==0))%聚群和追尾鱼太多太拥挤，都不发生；觅食觅不到更好的，造成三种行为都不发生。
海底地形测量图的插值模型MATLAB代码
%数据输入x_Yd=[129.0 140.0 108.5 88.0 185.5 195.5 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5];y_Yd=[ 7.5 141.5 28.0 147.0 22.5 137.5 85.5 -6.5 -81.0 3.0 56.5 -66.5 84.0 -38.5];z_Ft=[ 4 8 6 8 6 8 8 9 9 8 8 9 4 9];%由已知14点计算出14*14点（其中有两点重复）%计算QiQj的值Qi_QjQi_Qj=zeros(14,14);for i=1:14for j=1:14Qi_Qj(i,j)=sqrt((x_Yd(i)-x_Yd(j)).^2+(y_Yd(i)-y_Yd(j)).^2);endend[G_x,G_y]=meshgrid(x_Yd,y_Yd); %计算所求196个G点的坐标%求出每个G点对应的深度Z_gfor m=1:14for k=1:14if ((m==k)|((m==12)&&(k==13))|((m==13)&&(k==12))) %排除Qi,Qj为同一点和重复点的情况Z_g(m,k)=z_Ft(m);continue;else%计算GQij=min{GQi,GQj}G_Qi=zeros(1,14);G_Qj=zeros(1,14);for i=1:14G_Qi(i)=sqrt((G_y(m,k)-y_Yd(i)).^2+(G_x(m,k)-x_Yd(i)).^2);endfor j=1:14G_Qj(j)=sqrt((G_x(m,k)-x_Yd(j)).^2+(G_y(m,k)-y_Yd(j)).^2);endG_Qij=zeros(14,14);[G_Qi,G_Qj]=meshgrid(G_Qi,G_Qj);G_Qij=min(G_Qi,G_Qj);%计算G点到Qi,Qj两点连线的距离G_Pijfor i=1:14for j=1:14if sqrt((y_Yd(i)-y_Yd(j)).^2+(x_Yd(i)-x_Yd(j)).^2)G_Pij(i,j)=abs((G_x(m,k)-x_Yd(i)).*(y_Yd(i)-y_Yd(j))-(G_y(m,k)-y_Yd(i)).*(x_Yd(i)-x_Yd(j))). /sqrt((y_Yd(i)-y_Yd(j)).^2+(x_Yd(i)-x_Yd(j)).^2);endendend%计算PQij=min{PQi,PQj}P_Qi=sqrt(G_Qi.^2-G_Pij.^2);P_Qj=sqrt(G_Qj.^2-G_Pij.^2);P_Qij=min(P_Qi,P_Qj);% Z_g_wang=zeros(14,14);%根据平面几何关系计算出由Qi,Qj点的深度线性外推的Pij的深度Z_pfor i=1:14;for j=1:14;if (Qi_Qj(i,j)==0)|(z_Ft(i)==z_Ft(j))Z_p(i,j)=z_Ft(i);else if ((P_Qi(i,j)<=P_Qj(i,j))&&(Qi_Qj(i,j)<=P_Qi(i,j))&&(z_Ft(i)>z_Ft(j)))Z_p(i,j)=z_Ft(j)-(abs(z_Ft(i)-z_Ft(j))).*P_Qij(i,j)./Qi_Qj(i,j);else if ((P_Qi(i,j)<P_Qj(i,j))&&(Qi_Qj(i,j)<P_Qj(i,j))&&(z_Ft(j)>z_Ft(i)))Z_p(i,j)=abs(z_Ft(i)-(abs(z_Ft(j)-z_Ft(i)).*P_Qij(i,j)./Qi_Qj(i,j)));else if ((P_Qj(i,j)<=P_Qi(i,j))&&(Qi_Qj(i,j)<=P_Qi(i,j))&&(z_Ft(j)>z_Ft(i)))Z_p(i,j)=z_Ft(j)+(abs(z_Ft(j)-z_Ft(i)).*P_Qij(i,j)./Qi_Qj(i,j));else if ((P_Qi(i,j)<=P_Qj(i,j))&&(Qi_Qj(i,j)<=P_Qj(i,j))&&(z_Ft(i)>z_Ft(j)))Z_p(i,j)=z_Ft(i)+(abs(z_Ft(i)-z_Ft(j))).*P_Qij(i,j)./Qi_Qj(i,j);else if ((P_Qi(i,j)<=Qi_Qj(i,j))&&(P_Qj(i,j)<=Qi_Qj(i,j))&&(z_Ft(i)>z_Ft(j)))Z_p(i,j)=z_Ft(j)+(abs(z_Ft(i)-z_Ft(j))).*P_Qj(i,j)./Qi_Qj(i,j);else if ((P_Qi(i,j)<=Qi_Qj(i,j))&&(P_Qj(i,j)<=Qi_Qj(i,j))&&(z_Ft(j)>z_Ft(i)))Z_p(i,j)=z_Ft(i)+(abs(z_Ft(j)-z_Ft(i)).*P_Qi(i,j)./Qi_Qj(i,j));elseZ_P(i,j)=1;endendendendendendendendend%计算G点深度Z_gZ_g_1=ones(14,14)./(G_Pij.^2+G_Qij.^2+Qi_Qj.^2);Z_g_2=Z_p./(G_Pij.^2+G_Qij.^2+Qi_Qj.^2);for i=1:14for j=1:14Z_g_3(i,j)=sum(sum(Z_g_1));endendZ_g(m,k)=sum(sum(Z_g_2./Z_g_3));endendend%程序一：插值并作海底曲面图x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=-griddata(G_x,G_y,Z_g,x1,y1,'v4');meshc(x1,y1,z1)xlabel('x'),ylabel('y'),zlabel('z')figure;%程序二：插值并作出水深小于5的海域范围。
海杂波matlab仿真程序
海杂波MATLAB仿真程序介绍海洋中存在着各种各样的波浪，包括正弦波、涌浪、破碎波等。
其中，海杂波指的是由多个波浪相互作用形成的复杂波浪。
海杂波的研究对于海洋工程、海上交通以及海洋环境保护都具有重要的意义。
为了更好地理解和研究海杂波，我们可以利用MATLAB编写仿真程序，模拟出不同情况下的海杂波。
学习目标•了解海杂波的特点和形成原理•掌握使用MATLAB进行波浪模拟的方法和技巧•学会分析和处理仿真结果，得出相关结论和应用海杂波的特点和形成原理1. 特点海杂波具有以下特点： - 复杂性：海杂波是由多个波浪相互叠加形成的，其形状和波动规律都十分复杂。
- 非线性：海杂波的形成涉及到波浪的非线性相互作用，其振幅会发生剧烈变化。
- 随机性：海洋环境是一个复杂的系统，海杂波的形成和发展都具有一定的随机性。
2. 形成原理海杂波的形成原理主要包括以下几个方面： 1. 非线性作用：波浪在传播过程中会发生非线性相互作用，即波浪之间相互影响，使得波浪的振幅和形状发生变化。
这种非线性作用会导致波浪的增长和剧烈变化，进而形成海杂波。
2. 断裂和合并：在海洋中存在着各种不同频率和振幅的波浪。
当这些波浪相互碰撞时，会发生波浪的断裂和合并现象，使得波浪的形状发生变化，从而形成海杂波。
3. 激发机制：海洋中的风力、地震和潮汐等因素都会激发波浪的形成和发展。
这些外部力量的作用使得海杂波的形成更加复杂和多样化。
MATLAB中的海杂波仿真程序为了模拟海杂波的形成和发展过程，我们可以利用MATLAB编写仿真程序。
下面介绍一种基于MATLAB的海杂波仿真程序设计方法。
1. 数据初始化首先，我们需要初始化一些数据，包括海平面高度、时间步长、模拟时间等。
这些数据将作为仿真程序的输入。
2. 随机波浪生成在海洋中存在着各种不同频率和振幅的波浪。
我们可以利用MATLAB中的随机数生成函数来生成这些波浪，并模拟其在海洋中的传播和相互作用。
3. 波浪传播模拟利用MATLAB中的波动方程模型和数值求解方法，我们可以模拟波浪在海洋中的传播过程。
自然水系(或分水岭)自动提取的matlab实现方案
水系提取算法（D8）的matlab 实现D8算法是当今非常成熟的提取水系（或分水岭）的计算机实现方法。
这里笔者结合自身的编程实践介绍D8算法的matlab 实现过程，谨与大家分享。
D8算法分为两部分实现：每个栅格点水流方向的计算；每个栅格点汇水面积的计算。
1、水流方向的计算在水流方向的计算中，本人结合matlab 内嵌的滤波工具，利用一个简单的非线性滤波来实现，具体内容包括一个自己编辑的确定水流方向的函数flowdirection 和若干脚本命令，如下：以下是自己编辑的flowdirection 函数代码：function output=flowdirection(a)n=size(a);for i=1:n(2) k=0; b=-inf; for j=[1:4,6:9]if rem(j,2)==0r(j)=a(5,i)-a(j,i);elser(j)=(a(5,i)-a(j,i))/sqrt(2);endif r(j)>bb=r(j);k=j;endendoutput(i)=k;end执行的脚本代码如下：>> uiopen('F:\课件\2009田淑芳\正射校正\DEM.tif',1)%打开DEM 高程数据，生成的数据矩阵名可自行定义，这里默认为DEM>>DEM=mat2gray(DEM);%高程数据归一化处理>>DEM=padarray(DEM,[1 1],'replicate');%对数据边缘进行自动填充>>direc8=colfilt(DEM,[3 3],'sliding',@flowdirection);%执行非线性滤波，得到方向矩阵，得到的值分别为{1 ,2,3,4,6,7,8,9}.其方向意义如上图矩阵所示，如1代表北西向，8代表正东1、汇水面积的计算这里笔者根据自己的编程实践经验，给出一个函数来计算不同窗口大小下汇水面积的大小。
生成二维pm谱海洋面的matlab程序代码
生成二维PM谱海洋面的Matlab程序代码1. 介绍PM(Phased-Mission)谱海洋是指在多任务执行过程中的频谱特性，它在通信、雷达、遥感等领域有着广泛的应用。
生成二维PM谱海洋面的Matlab程序代码可以帮助工程师们分析并优化系统的性能，因此具有重要的研究和应用价值。
2. PM谱海洋的定义PM谱海洋是指在多任务执行过程中，系统任务性能不均匀变化所呈现出的频谱特性。
可以描述为任务完整度随时间和频率的变化。
在实际应用中，PM谱海洋通常为二维图像，横轴代表频率，纵轴代表时间。
3. 生成二维PM谱海洋面的Matlab程序代码接下来我们将介绍生成二维PM谱海洋面的Matlab程序代码，以供工程师们学习和参考。
以下是完整的Matlab程序代码：```matlab生成二维PM谱海洋面的Matlab程序代码参数设置fs = 1000; 采样频率T = 10; 采样时长t = 0:1/fs:T; 时间序列f0 = 100; 信号频率f1 = 200; 信号频率增量N = length(t); 信号长度生成信号x = sin(2*pi*f0*t) + sin(2*pi*(f0+f1)*t);信号频谱X = fft(x,N);f = fs/N*(0:N-1);P = abs(X).^2/N;生成PM谱海洋[pxx,f,t] =pspectrum(x,fs,'spectrogram','FrequencyLimits',[0,400]); mesh(t,f,10*log10(pxx)); 画出PM谱海洋面xlabel('Time(s)'); 横轴标签ylabel('Frequency(Hz)'); 纵轴标签zlabel('Power/Frequency(dB/Hz)'); z轴标签title('PM Spectrogram'); 图像标题```4. 使用说明以上代码通过调用Matlab中的pspectrum函数，可以生成二维PM谱海洋面。
matlab 风驱动 海水垂向混合计算
matlab 风驱动海水垂向混合计算摘要：一、引言1.风驱动海水混合的重要性2.海水垂向混合计算的研究背景二、Matlab 简介1.Matlab 的基本功能2.Matlab 在风驱动海水混合计算中的应用三、海水垂向混合计算方法1.垂向混合模型2.数值计算方法四、Matlab 实现海水垂向混合计算1.编写计算程序2.输入参数设置3.运行结果分析五、计算结果与分析1.风驱动海水混合的数值模拟2.结果验证与讨论六、结论1.Matlab 在风驱动海水混合计算中的优势2.对未来研究的展望正文：一、引言风驱动海水混合在海洋环境中具有重要作用，它影响着海洋生态系统的稳定性、养分循环以及气候模式等。
随着科技的发展，人们越来越关注如何更准确地模拟和预测风驱动海水混合过程。
Matlab 作为一种功能强大的数学软件，已被广泛应用于各种科学计算领域。
本文将介绍如何利用Matlab 进行风驱动海水垂向混合计算，以期为相关研究提供参考。
二、Matlab 简介Matlab 是一种商业数学软件，由美国MathWorks 公司开发。
它具有丰富的数学函数库、图形绘制功能以及强大的编程能力，可以进行各种科学计算、数据分析、建模以及可视化等。
在风驱动海水混合计算中，Matlab 可以用于建立数值模型、求解偏微分方程以及分析计算结果等。
三、海水垂向混合计算方法为了模拟风驱动海水混合过程，首先需要建立垂向混合模型。
一般采用Navier-Stokes 方程和能量守恒方程描述海水混合过程。
在数值计算方面，通常采用有限差分法、有限元法等方法对偏微分方程进行离散求解。
四、Matlab 实现海水垂向混合计算本文以一个简单的风驱动海水垂向混合模型为例，介绍如何利用Matlab 进行计算。
首先，根据建立的垂向混合模型，编写相应的计算程序。
接着，设置输入参数，如风速、海水密度等。
最后，运行程序并获得计算结果。
通过对结果的分析，可以了解风驱动海水混合过程中的动力学和热力学特征。
人工水母搜索算法—matlab代码
⼈⼯⽔母搜索算法—matlab代码clcclearfoj = @ Sphere;Lb = -100; % 搜索空间下界Ub = 100; % 搜索空间上界N_iter = 1000; % 最⼤迭代次数n_pop = 50; % 种群个数d = 10; % 种群维度beta = 3;gamma = 0.1;Z = zeros(n_pop, d);% 随机⽣成⼀个d维向量Z(1, :) = rand(1, d);% 利⽤logistic⽣成n_pop个向量for i=2:n_popZ(i,:) = 4.0*Z(i-1,:).*(1-Z(i-1,:));end% 将z的各个分量载波到对应变量的取值区间pop = zeros(n_pop, d);fitness = zeros(n_pop, 1);for i=1:n_poppop(i,:) = Lb + (Ub - Lb)*Z(i,:);fitness(i) = foj(pop(i,:));endmu = mean(pop, 1);figurescatter(pop(:,1), pop(:,2), 'r*');hold onscatter(mu(:,1), mu(:,2),'ko');box on[bestScore, loc] = min(fitness);bestPop = pop(loc, :); % 当前种群中⾷物数⽬最多的位置S = pop; % S是更新后的种群fmin = zeros(N_iter,1);for t=1:N_iterfor i=1:n_popc = abs((1-t/N_iter)*(2*rand-1));if c>=0.5% 洋流位置更新公式trend = bestPop - beta*rand*mu; % 洋流的⽅向,即论⽂中的公式.9S(i,:) = pop(i,:) + rand(1,d).*trend; % ⽔母位置更新，即论⽂中的公式.11else% 种群内部运动if rand>(1-c)S(i,:) = pop(i,:) + gamma*rand(1,d).*(Ub - Lb); % 被动运动，即论⽂公式.12elseJK = randperm(n_pop);if foj(pop(JK(1)))>=foj(pop(JK(2)))Direction = pop(JK(2),:) - pop(JK(1),:);elseDirection = pop(JK(1),:) - pop(JK(2),:);endStep = rand(1, d).*Direction;S(i,:) = pop(i,:) + Step; % 主动运动，论⽂公式.16endend% 检查边界S(i,:) = simpleBound(S(i,:), Lb, Ub);Fnew = foj(S(i,:));% 判断是否更新相应位置和适应度值if Fnew<fitness(i)fitness(i) = Fnew;pop(i,:) = S(i,:);end% 判断是否更新最优值if Fnew<bestScorebestScore = Fnew;bestPop = S(i,:);endendfmin(t) = bestScore;% 命令窗⼝输出disp(['Iteration ' num2str(t) ': Best Cost = ' num2str(bestScore)]);% disp(['Bestpop ' num2str(bestPop)]);endfiguresemilogy(fmin, 'r-.')function x=simpleBound(x,Lb,Ub) % 函数名称：越界处理函数 % param x：⽔母 % param Lb：变量下界 % param Ub：变量上界 d = length(x); for i=1:d if x(i) > Ub x(i) = (x(d) - Ub) + Lb; elseif x(i) < Lb x(i) = (x(d) - Lb) + Ub; end endendfunction [y] = Sphere(xx)d = length(xx);sum = 0;for ii = 1:dxi = xx(ii);sum = sum + xi^2;endy = sum;end。
用Matlab软件计算
用Matlab 软件编程计算1、分别用逆矩阵法、LU 分解和QR 分解解线性方程组：⎪⎪⎩⎪⎪⎨⎧-=+++=+++-=+++=---2230245132744321432143214321x x x x x x x x x x x x x x x x 。
2、分别用求特征值法和直接法求方程0432324=+--x x x 的全部解。
3、输入a 、b 、c 的值，求一元二次方程)0(02≠=++a c bx x a 的解。
4、某商场对顾客所购买的商品实行打折销售，标准如下(商品价格用price 来表示)：price<100 没有折扣 100≤price<400 3%折扣 400≤price<1000 5%折扣 1000≤price<2000 8%折扣 2000≤price<4000 10%折扣4000≤price 14%折扣试设计一个程序，输入所售商品的价格，输出实际销售价格。
5、求∑=+1001121i i 的结果。
6、利用函数文件，实现将极坐标(ρ,θ)（θ单位：度，deg ）转换成直角坐标(x,y)。
7、Fibonacci 数列定义如下：⎪⎩⎪⎨⎧>+===--)2(112121n ff f f f n n n，试用递归调用函数求Fibonacci 数列的第10项。
8、分别以条形图、阶梯图、杆图和填充图形式绘制曲线y=xcos(2x)在π≤≤x 0上的图象。
9、利用子图作出隐函数01622=-+y x ，0422=--y x ，04=-+xy xe y 的图象。
10、在区间]4,0[π上画出参数曲线x=sin(2t)，y=cos(2t)，z=t ，并分别标注坐标轴。
11、对矩阵Ｘ＝⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-----324.325.235.192.021.35.223.298.101.29.32.33.15.2，试求访矩阵的特征值和特征向量，并从不同维方向求平均值和标准方差。
海水密度计算matlab
海水密度计算matlab
在MATLAB中，可以使用以下公式来计算海水的密度：
海水密度 = 海水的标准密度 [1 (海水的温度 4) 0.0005 + (海水的盐度 35) 0.001]
其中，海水的标准密度为 1025 kg/m^3，海水的温度以摄氏度
为单位，盐度以盐度计为单位。
在MATLAB中，你可以编写一个函数来计算海水的密度，例如：
matlab.
function density = calculateSeaWaterDensity(temperature, salinity)。
standardDensity = 1025;
density = standardDensity (1 (temperature 4)
0.0005 + (salinity 35) 0.001);
end.
然后你可以调用这个函数，并传入海水的温度和盐度来计算海水的密度，例如：
matlab.
temperature = 10; % 海水的温度为10摄氏度。
salinity = 36; % 海水的盐度为36盐度计。
density = calculateSeaWaterDensity(temperature, salinity);
disp(['海水的密度为，', num2str(density), ' kg/m^3']);
这样就可以在MATLAB中计算海水的密度了。
当然，这只是一个简单的示例，实际应用中可能还需要考虑更多因素，比如压力等。
希望这个回答能够帮到你。
第4讲matlab海洋大气数据可视化
例：whos plot(x1, y1, x2, y2); figure;plot(x1,y,x,y)
3．具有两个纵坐标标度的图形 在MATLAB中，如果需要绘制出具有不同纵坐标标 度的两个图形，可以使用plotyy绘图函数。调用 格式为： plotyy(x1,y1,x2,y2) 其中x1,y1对应一条曲线，x2,y2对应另一条曲线。横 坐标的标度相同，纵坐标有两个，左纵坐标用于 x1,y1数据对，右纵坐标用于x2,y2数据对。
doc; latex letter
例5-7 在0≤x≤2区间内，绘制曲线y1=2e-0.5x和 y2=cos(4πx)，并给图形添加图形标注。 程序如下： x=0:pi/100:2*pi; y1=2*exp(-0.5*x); y2=cos(4*pi*x); plot(x,y1,x,y2) title('x from 0 to 2{\pi}'); %加图形标题 xlabel('Variable X'); %加X轴说明 ylabel('Variable Y'); %加Y轴说明 text(0.8,1.5,'曲线y1=2e^{-0.5x}'); %在指定位置 添加图形说明 text(2.5,1.1,'曲线y2=cos(4{\pi}x)'); legend(‘y1’,‘ y2’) %加图例
5.1.6 图形窗口的分割 subplot函数的调用格式为： subplot(m,n,p) 该函数将当前图形窗口分成m×n个绘图区， 即每行n个，共m行，区号按行优先编号， 且选定第p个区为当前活动区。在每一个绘 图区允许以不同的坐标系单独绘制图形。 例5-10 在图形窗口中，以子图形式同时绘制 多根曲线。
5.1.7 二维统计分析图 在MATLAB中，二维统计分析图形很多，常 见的有条形图、阶梯图、杆图和填充图等， 所采用的函数分别是： bar(x,y,选项) stairs(x,y,选项) stem(x,y,选项) fill(x1,y1,选项1,x2,y2,选项2,…)
matlab 海温区域平均
matlab 海温区域平均
【最新版】
目录
1.Matlab 简介
2.海温区域平均的概念
3.Matlab 在海温区域平均计算中的应用
4.具体操作步骤
5.总结
正文
【1.Matlab 简介】
Matlab 是一种广泛应用于科学计算和数据分析的编程语言。
它以其
高效的矩阵计算和强大的可视化功能，成为科研和工程领域中的重要工具。
【2.海温区域平均的概念】
海温区域平均是指在一定地理范围内，海洋表面温度的平均值。
这个值可以反映出该区域的海洋温度状况，是研究海洋环流、海洋生态和气候变化的重要参数。
【3.Matlab 在海温区域平均计算中的应用】
Matlab 可以用来方便地计算海温的区域平均。
用户只需输入几行代码，就可以完成大量的数据处理和计算任务。
【4.具体操作步骤】
以下是使用 Matlab 计算海温区域平均的具体操作步骤：
1) 首先，需要安装并打开 Matlab 软件。
2) 然后，读取海温数据。
这些数据可以是来自于网络或者其他来源
的。
3) 使用 Matlab 的区域函数，对海温数据进行分区域平均计算。
4) 最后，将计算结果进行可视化，以便于观察和分析。
【5.总结】
总的来说，Matlab 在海温区域平均计算中的应用，极大地提高了数据处理的效率和精确度。
matlab 淹膜海洋函数
matlab 淹膜海洋函数
Matlab淹膜海洋函数是一种用于模拟海洋环境的工具，它能够模拟海洋中的各种物理现象，包括波浪、水流、海底地形等等。
通过使用该函数，可以进行海洋工程、海洋资源开发和海洋环境改善等方面的研究。
Matlab淹膜海洋函数的使用非常简单，只需要输入一些基本的参数，例如海洋深度、波浪高度、流速等等，即可生成海洋环境的模拟结果。
此外，Matlab淹膜海洋函数还提供了一些可视化工具，可以帮助用户更直观地分析海洋环境的各种特征。
总之，Matlab淹膜海洋函数是一种非常强大的工具，可以帮助用户更好地理解海洋环境的特点，为海洋工程和海洋资源开发提供有力的支持。
- 1 -。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
代码：clear all;close all;clc%M2,S2,N2,K2,K1,O1,P1,Q1八个分潮的杜德深常数U1、2、3、4、5、6、0U=[2,0,0,0,0,0,0;2,2,-2,0,0,0,0;2,-1,0,1,0,0,0;2,2,0,0,0,0,0;1,1,0,0,0,0,1;1,-1,0,0,0,0,-1;1,1,-2,0,0,0,-1;1,-2,0,1,0,0,-1;];%% 天文元素的变化速度twys_[]twys_=[14.4921,0.549,0.0411,0.0046,0.0022,0.0000000196].*(pi/180); %% 求天文要素twys[tao,s,h_,p,N_,p_]year=2003;mon=3;day=3+15;time=0;i=floor((year-1900)/4);ddd=76;dy=year-1900;di=ddd+i+time/24;s=277.02+129.3848*dy+13.1764*di;h_=280.19-0.2387*dy+0.9857*di;p=334.39+40.6625*dy+0.1114*di;N_=100.84+19.2382*dy+0.053*di;p_=281.22+0.0172*dy+0.00005*di;tao=15*time-s+h_;twys=[tao,s,h_,p,N_,p_].*pi/180;%% 求[M2,S2,N2,K2,K1,O1,P1,Q1]八个分潮角速度sig、初始相位V sigg=zeros(1,8);V0=zeros(1,8);for k =1:8sigg(k)=sum(U(k,1:6).*twys_);V0(k)=sum(U(k,1:7).*[twys,pi/2]);V0(k)=rem(V0(k),(2*pi));if abs(V0(k))<0.0001V0(k)=0;endif V0(k)<0V0(k)=V0(k)+2*pi;endend%% 求交点因子ff、交点订正角uufK1_cosuK1=1+0.002*cos(-2*twys(4)-twys(5))+0.0001*cos(-2*twys(5)) -0.0198*cos(-twys(5))+...0.1356*cos(twys(5))-0.0029*cos(2*twys(5));fK1_sinuK1=0.002*sin(-2*twys(4)-twys(5))+0.0001*sin(-2*twys(5))-0.0 198*sin(-twys(5))+...0.1356*sin(twys(5))-0.0029*sin(2*twys(5));uK1=atan(fK1_sinuK1/fK1_cosuK1)+2*pi;fK1=abs(sqrt(fK1_cosuK1^2+fK1_sinuK1^2));%%%计算M2分潮的交点因子和交点订正角fM2_cosuM2=1+0.0005*cos(-2*twys(5))-0.0373*cos(-twys(5))+0.0006* cos(2*twys(4))+...0.0002*cos(2*twys(4)+twys(5));fM2_sinuM2=0.0005*sin(-2*twys(5))-0.0373*sin(-twys(5))+0.0006*sin (2*twys(4))+...0.0002*sin(2*twys(4)+twys(5));uM2=atan(fM2_sinuM2/fM2_cosuM2)+2*pi;fM2=abs(sqrt(fM2_cosuM2^2+fM2_sinuM2^2));%计算K2分潮的交点因子和交点订正角fK2_cosuK2=1-0.0128*cos(-twys(5))+0.2890*cos(twys(5))+0.0324*cos (2*twys(5));fK2_sinuK2=-0.0128*sin(-twys(5))+0.2890*sin(twys(5))+0.0324*sin(2*t wys(5));uK2=atan(fK2_sinuK2/fK2_cosuK2)+2*pi;fK2=abs(sqrt(fK2_cosuK2^2+fK2_sinuK2^2));%%%计算P1分潮的交点因子和交点订正角fP1_cosuP1=1+0.0008*cos(-2*twys(5))-0.0112*cos(-twys(5))-0.0015*co s(2*twys(4))...-0.0003*cos(2*twys(4)+twys(5));fP1_sinuP1=0.0008*cos(-2*twys(5))-0.0112*cos(-twys(5))-0.0015*cos(2 *twys(4))...-0.0003*cos(2*twys(4)+twys(5));uP1=atan(fP1_sinuP1/fP1_cosuP1)+2*pi;fP1=abs(sqrt(fP1_cosuP1^2+fP1_sinuP1^2));%%%计算O1分潮的交点因子和交点订正角fO1_cosuO1=1-0.0058*cos(-2*twys(5))+0.1885*cos(-twys(5))+0.0002*c os(2*twys(4)-twys(5))...-0.0064*cos(2*twys(4))-0.0010*cos(2*twys(4)+twys(5));fO1_sinuO1=-0.0058*sin(-2*twys(5))+0.1885*sin(-twys(5))+0.0002*sin (2*twys(4)-twys(5))...-0.0064*sin(2*twys(4))-0.0010*sin(2*twys(4)+twys(5));uO1=atan(fO1_sinuO1/fO1_cosuO1);fO1=abs(sqrt(fO1_cosuO1^2+fO1_sinuO1^2));fS2=1;uS2=0;fQ1=fO1;uQ1=uO1;fN2=fM2;uN2=uM2;ff=[fM2,fS2,fN2,fK2,fK1,fO1,fP1,fQ1];ff=[0.9837,1,0.9837,1.511,1.0632,0.0968,0.9936,1.0968];uu=[uM2,uS2,uN2,uK2,uK1,uO1,uP1,uQ1];uu=[6.2504,0,6.2504,6.0167,6.1549,0.144,6.2724,0.144];%% 最小二乘法计算法方程date=load('E:\magicalcity\2020.3\海洋要素计算\上机实验\3.txt'); %构造系数矩阵AN=721;N_=360;A=zeros(9,9);for f=2:9aa=(sin(sigg(f-1)*N/2));bb=(sin(sigg(f-1)/2));A(1,f)=aa/bb;A(f,1)=A(1,f);endA(1,1)=N;for k=2:9for j=2:9if k==jA(k,j)=(N+sin(sigg(k-1)*N)/sin(sigg(k-1)))/2;endif k~=ja=(sin((sigg(k-1)-sigg(j-1))*N/2));b=(sin((sigg(k-1)-sigg(j-1))/2));c=(sin((sigg(k-1)+sigg(j-1))*N/2));d=(sin((sigg(k-1)+sigg(j-1))/2));A(k,j)=(a/b+c/d)/2;endendend%构造系数矩阵BB=zeros(8,8);for k=1:8for j=1:8if k==jB(k,j)=(N-sin(sigg(k)*N)/sin(sigg(k)))/2;endif k~=ja=(sin((sigg(k)-sigg(j))*N/2));b=(sin((sigg(k)-sigg(j))/2));c=(sin((sigg(k)+sigg(j))*N/2));d=(sin((sigg(k)+sigg(j))/2));B(k,j)=(a/b-c/d)/2;endendend%构造F1、F2F1=zeros(9,1);F2=zeros(8,1);F1(1,1)=sum(date(:,2));for k=2:9for j=1:NF1(k,1)=F1(k,1)+date(j,2)*cos(sigg(k-1)*(j-361));F2(k-1,1)=F2(k-1,1)+date(j,2)*sin(sigg(k-1)*(j-361));endend%法方程的解x，yx=A\F1;y=B\F2;%% 计算准调和振幅R和位相sitafor k=2:9R(k-1)=(x(k)*x(k)+y(k-1)*y(k-1))^(1/2);if x(k)>=0sita(k-1)=asin(y(k-1)/R(k-1));endif x(k)<0sita(k-1)=pi-asin(y(k-1)/R(k-1));endend%% 计算调和常数H、gfor k=1:8H(k)=R(k)./ff(k);g(k)=sita(k)+V0(k)+uu(k);g(k)=rem(g(k),2*pi);if g(k)<0g(k)=g(k)+2*pi;endg(k)=g(k)*180/pi;end%% 潮汐自报hhdelt=0;for k=1:Nhh(k)=0;for j=1:8hh(k)=hh(k)+ff(j)*H(j)*cos(sigg(j)*(k-361)+V0(j)+uu(j)-g(j));endhh(k)=hh(k)+x(1);delt=delt+(hh(k)-date(k,2))^2;enddelt=(delt/N)^(1/2);figure(1)plot(date(:,2),'r')hold on;plot(hh(1,:),'b')title('实测（红色）、自报（蓝色）潮汐曲线','fontsize',12)%% 5.1-6.1的潮汐预报hhpfor k=1:745hhp(k)=0;for j=1:8hhp(k)=hhp(k)+ff(j)*H(j)*cos(sigg(j)*(k+1056)+V0(j)+uu(j)-g (j));endendfigure(2)plot(hhp(1,:),'b')title('预报5.1-6.1潮汐','fontsize',12)。