%
%--------------------------------------------------------------------------------------------!
% 1D CSAMT
%
%--------------------------------------------------------------------------------------------!
function Source1(varargin)
clear all; %clear functions;
global unit2fid; if ~isempty(unit2fid), unit2fid=[]; end
persistent ex ey ez frq frqname frquencecount hx hy hz i modelname mu0 pi r rou sourcecoordinatename surveycoordinatename surveycount x0 x1 y0 y1 ;
if isempty(surveycount), surveycount=1; end;
if isempty(frquencecount), frquencecount=41; end;
if isempty(pi), pi=3.141592654; end;
if isempty(mu0), mu0=4.0.*pi.*1.0e-7; end;
if isempty(i), i=0; end;
if isempty(x0), x0=zeros(1,2); end;
if isempty(y0), y0=zeros(1,2); end;
if isempty(x1), x1=zeros(1,surveycount); end;
if isempty(y1), y1=zeros(1,surveycount); end;
if isempty(rou), rou=0; end;
if isempty(r), r=zeros(1,surveycount); end;
if isempty(frq), frq=zeros(1,frquencecount); end;
if isempty(sourcecoordinatename), sourcecoordinatename=repmat(' ',1,20); end;
if isempty(surveycoordinatename), surveycoordinatename=repmat(' ',1,20); end;
if isempty(modelname), modelname=repmat(' ',1,20); end;
if isempty(frqname), frqname=repmat(' ',1,20); end;
if isempty(ex), ex=zeros(1,frquencecount); end;
if isempty(ey), ey=zeros(1,frquencecount); end;
if isempty(ez), ez=zeros(1,frquencecount); end;
if isempty(hx), hx=zeros(1,frquencecount); end;
if isempty(hy), hy=zeros(1,frquencecount); end;
if isempty(hz), hz=zeros(1,frquencecount); end;
[writeErrFlag]=writeFmt(1,'%c','''please input SourceCoordinate''');
sourcecoordinatename=strAssign(sourcecoordinatename,[],[],input('','s'));
thismlfid=fopen(strtrim(sourcecoordinatename),'r+');
unit2fid=[unit2fid;11,thismlfid];
for i=1:2;
[readErrFlag,readEndFlag]=readFmt(11,['%g','%g'],'x0(i)','y0(i)');
end; %i=fix(2+1);
[writeErrFlag]=writeFmt(1,['%c'],'''please input SurveyCoordinate''');
surveycoordinatename=strAssign(surveycoordinatename,[],[],input('','s'));
thismlfid=fopen(strtrim(surveycoordinatename),'r+');
unit2fid=[unit2fid;12,thismlfid];
for i=1:surveycount;
[readErrFlag,readEndFlag]=readFmt(12,['%g','%g'],'x1(i)','y1(i)');
end; i=fix(surveycount+1);
[writeErrFlag]=writeFmt(1,['%c'],'''please input frq information''');
frqname=strAssign(frqname,[],[],input('','s'));
thismlfid=fopen(strtrim(frqname),'r+');
unit2fid=[unit2fid;13,thismlfid];
for i=1:frquencecount;
[readErrFlag,readEndFlag]=readFmt(13,['%g'],'frq(i)');
end; i=fix(frquencecount+1);
[writeErrFlag]=writeFmt(1,['%c'],'''please input model''');
modelname=strAssign(modelname,[],[],input('','s'));
thismlfid=fopen(strtrim(modelname),'r+');
unit2fid=[unit2fid;14,thismlfid];
[readErrFlag,readEndFlag]=readFmt(14,['%g'],'rou');
%计算收发距
for i=1:surveycount;
r(i)=sqrt(((x0(2)-x0(1))./2.0-x1(i)).^2+((y0(2)-y0(1))./2.0-y1(i)).^2);
end; i=fix(surveycount+1);
%对每一个测点都计算一次场值
for i=1:surveycount;
[frquencecount,mu0,pi,x1(i),y1(i),r(i),rou,frq,ex,ey,ez,hx,hy,hz]=calculatefield(frquencecount,mu0,pi,x1(i),y1(i),r(i),rou,frq,ex,ey,ez,hx,hy,hz);
end; i=fix(surveycount+1);
thismlfid=fopen(strtrim('test.dat'),'r+');
unit2fid=[unit2fid;20,thismlfid];
for i=1:frquencecount;
%write(*,*) frq(i),(real(ex(i))**2+imag(ex(i))**2)/(real(hy(i))**2+imag(hy(i))**2)/5.0/frq(i)
[writeErrFlag]=writeFmt(20,['%g','%*g'],'frq(i)','real(ex(i))');
end; i=fix(frquencecount+1);
pause;
error(['Even though Matlab just threw an error, there may or may not be a problem with this matlab code. This message merely signifies that a stop was encountered at this location in the original fortran code from which this matlab code was transalted. Please check the code to decide whether this is an actual error condition, or (which is often the case) whether this "stop" in fortran merely indicates the end of natural execution of the code. ',char(10),';']);
end %program
%---------------------------------------------------------------------------------------------------------!
%FrquenceCount:频点个数 mu0:磁导率 siteX:测点横坐标 siteY:测点综坐标 r:收发距 frq:频率 rou:电阻率
%a0,a1,s0,s1: hankel滤波系数的参数 y0,y1:导纳 z0:阻抗 ii:电流强度 dl:导线长度 delta:电导率 vacuum:介电常数
%---------------------------------------------------------------------------------------------------------!
function [frquencecount,mu0,pi,sitex,sitey,r,rou,frq,ex,ey,ez,hx,hy,hz]=calculatefield(frquencecount,mu0,pi,sitex,sitey,r,rou,frq,ex,ey,ez,hx,hy,hz,varargin);
persistent a0 a1 delta dl eii ex1 ex2 ey1 ey2 ez1 ez2 firstCall hank0_120 hank1_140 hx1 hx2 hy1 hy2 hz1 hz2 i ii j m0 m1 rte rtm s0 s1 u0 u1 vacuum y0 y1 z0 ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(a0), a0=0; end;
if isempty(a1), a1=0; end;
if isempty(s0), s0=0; end;
if isempty(s1), s1=0; end;
if isempty(m0), m0=zeros(1,120); end;
if isempty(m1), m1=zeros(1,140); end;
if isempty(ex1), ex1=zeros(1,frquencecount); end;
if isempty(ex2), ex2=zeros(1,frquencecount); end;
if isempty(ey1), ey1=zeros(1,frquencecount); end;
if isempty(ey2), ey2=zeros(1,frquencecount); end;
if isempty(ez1), ez1=zeros(1,frquencecount); end;
if isempty(ez2), ez2=zeros(1,frquencecount); end;
if isempty(hx1), hx1=zeros(1,frquencecount); end;
if isempty(hx2), hx2=zeros(1,frquencecount); end;
if isempty(hy1), hy1=zeros(1,frquencecount); end;
if isempty(hy2), hy2=zeros(1,frquencecount); end;
if isempty(hz1), hz1=zeros(1,frquencecount); end;
if isempty(hz2), hz2=zeros(1,frquencecount); end;
if isempty(y0), y0=0; end;
if isempty(z0), z0=0; end;
if isempty(u0), u0=0; end;
if isempty(u1), u1=0; end;
if isempty(y1), y1=0; end;
if isempty(eii), eii=complex(0,1.0); end;
if isempty(rte), rte=0; end;
if isempty(rtm), rtm=0; end;
if isempty(ii), ii=1.0; end;
if isempty(dl), dl=1.0; end;
if isempty(delta), delta=0; end;
if isempty(vacuum), vacuum=8.854187817.*1.0e-12; end;
a0=-8.38850000000d+00;
s0=9.04226468670d-02;
a1=-7.91001919000d+00;
s1=8.79671439570d-02;
delta=1.0./rou;
if firstCall, [ hank0_120(1:120)]=[9.62801364263d-07,-5.02069203805d-06,1.25268783953d-05,-1.99324417376d-05,2.29149033546d-05,-2.04737583809d-05,1.49952002937d-05,-9.37502840980d-06,5.20156955323d-06,-2.62939890538d-06,1.26550848081d-06,-5.73156151923d-07,2.76281274155d-07,-1.09963734387d-07,7.38038330280d-08,-9.31614600001d-09,3.87247135578d-08,2.10303178461d-08,4.10556513877d-08,4.13077946246d-08,5.68828741789d-08,6.59543638130d-08,8.40811858728d-08,1.01532550003d-07,1.26437360082d-07,1.54733678097d-07,1.91218582499d-07,2.35008851918d-07,2.89750329490d-07,3.56550504341d-07,4.39299297826d-07,5.40794544880d-07,6.66136379541d-07,8.20175040653d-07,1.01015545059d-06,1.24384500153d-06,1.53187399787d-06,1.88633707689d-06,2.32307100992d-06,2.86067883258d-06,3.52293208580d-06,4.33827546442d-06,5.34253613351d-06,6.57906223200d-06,8.10198829111d-06,9.97723263578d-06,1.22867312381d-05,1.51305855976d-05,1.86329431672d-05,2.29456891669d-05,2.82570465155d-05,3.47973610445d-05,4.28521099371d-05,5.27705217882d-05,6.49856943660d-05,8.00269662180d-05,9.85515408752d-05,1.21361571831d-04,1.49454562334d-04,1.84045784500d-04,2.26649641428d-04,2.79106748890d-04,3.43716968725d-04,4.23267056591d-04,5.21251001943d-04,6.41886194381d-04,7.90483105615d-04,9.73420647376d-04,1.19877439042d-03,1.47618560844d-03,1.81794224454d-03,2.23860214971d-03,2.75687537633d-03,3.39471308297d-03,4.18062141752d-03,5.14762977308d-03,6.33918155348d-03,7.80480111772d-03,9.61064602702d-03,1.18304971234d-02,1.45647517743d-02,1.79219149417d-02,2.20527911163d-02,2.71124775541d-02,3.33214363101d-02,4.08864842127d-02,5.01074356716d-02,6.12084049407d-02,7.45146949048d-02,9.00780900611d-02,1.07940155413d-01,1.27267746478d-01,1.46676027814d-01,1.62254276550d-01,1.68045766353d-01,1.52383204788d-01,1.01214136498d-01,-2.44389126667d-03,-1.54078468398d-01,-3.03214415655d-01,-2.97674373379d-01,7.93541259524d-03,4.26273267393d-01,1.00032384844d-01,-4.94117404043d-
评论1