From 882604f974523fbde7291c94ac510499749de039 Mon Sep 17 00:00:00 2001
From: "andreas.lagg" <lagg@mps.mpg.de>
Date: Sat, 26 Oct 2019 09:22:00 +0200
Subject: [PATCH] added the support for using the atm_ outpout files as
ATM_INPUT files (requires the change of order in the fits file)
---
configure.ac | 2 +-
fortran/src/read_atminput.f90 | 31 ++++++++++++++++++++++++-------
fortran/src/read_fits4d.f90 | 2 +-
idlpro/atminput_fill.pro | 7 ++++---
idlpro/helix.pro | 1 +
idlpro/read_atminput.pro | 9 ++++++++-
6 files changed, 39 insertions(+), 13 deletions(-)
diff --git a/configure.ac b/configure.ac
index 805215b..f2cfa6f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -200,7 +200,7 @@ if test "x$with_mpi" == "x" ; then
[AC_MSG_ERROR([dislin.f90 or dislin.mod not found.])])])
HELIBS=$HELIBS" -ldislin"
if test `uname` == "Darwin" ; then
- HELIBS=$HELIBS" -lXm"
+ HELIBS=$HELIBS" -lXm -L/usr/local/lib/ "
fi
PPFLAGS=$PPFLAGS" "$PP"X11"
AC_MSG_NOTICE([X11 (DISLIN) version])
diff --git a/fortran/src/read_atminput.f90 b/fortran/src/read_atminput.f90
index f0ed9a0..de1c140 100644
--- a/fortran/src/read_atminput.f90
+++ b/fortran/src/read_atminput.f90
@@ -68,8 +68,8 @@ contains
case('ALPHA')
ipt%atm(ia)%par%ff=atminput(ia)%data(i,xv,yv)
case default
+ if (ipt%verbose.ge.1) &
write(*,*) 'Could not find parameter from ATMINPUT: ',atminput(ia)%name(i)
- call stop
end select
end if
end do
@@ -81,10 +81,12 @@ contains
use ipt_decl
use atminput_comm
integer(2) :: ia,i,nparmax=0
- integer(4) :: unit,npar,nx,ny,bsz,status,naxes(3),nfound
- integer(4) :: fpixels(3),lpixels(3),inc(3)
+ integer(4) :: unit,npar,nx,ny,bsz,status,naxes(3),fnaxes(3),nfound
+ integer(4) :: fpixels(3),lpixels(3),inc(3),fstatus
logical :: anyf
- ! real(4),allocatable :: data(:,:,:)
+ real(4),allocatable :: data(:,:,:)
+ character(LEN=72) :: comment
+ character(LEN=maxstr) :: code
allocate(atminput(ipt%ncomp))
@@ -102,7 +104,14 @@ contains
call stop
end if
!size & offset of old data
- call ftgknj(unit,'NAXIS',1,3,naxes,nfound,status)
+ call ftgknj(unit,'NAXIS',1,3,fnaxes,nfound,status)
+ naxes=fnaxes
+ !check if file is atm_ file from HeLIx+ and transpose data accordingly
+ call ftgkys(unit,'CODE',code,comment,fstatus)
+ if (trim(code).eq.'HELIX+') then
+ naxes(1)=fnaxes(3)
+ naxes(2:3)=fnaxes(1:2)
+ end if
if (ipt%verbose.ge.2) write(*,*) 'npar=',naxes(1),' nx=',naxes(2),' ny=',naxes(3)
if (naxes(1).gt.nparmax) nparmax=naxes(1)
if (ia.eq.1) then
@@ -121,17 +130,25 @@ contains
atminput(ia)%nx=nx
atminput(ia)%ny=ny
allocate(atminput(ia)%data(npar,nx,ny))
+ allocate(data(fnaxes(1),fnaxes(2),fnaxes(3)))
status=0
call ftgkns(unit,'PAR',1,npar,atminput(ia)%name(1:npar),nfound,status)
status=0
fpixels(:)=(/1,1,1/)
inc(:)=1
- lpixels=naxes
- call ftgsve(unit,1,3,naxes,fpixels,lpixels,inc,0,atminput(ia)%data,anyf,status)
+ lpixels=fnaxes
+ call ftgsve(unit,1,3,fnaxes,fpixels,lpixels,inc,0,data,anyf,status)
if (status.ne.0) then
write(*,*) 'Could not read data from ATM_INPUT file.'
call stop
end if
+ if (trim(code).eq.'HELIX+') then
+ atminput(ia)%data=reshape(data,shape(atminput(ia)%data),order=[2,3,1])
+ else
+ atminput(ia)%data=data
+ end if
+
+ deallocate(data)
call ftclos(unit,status)
end if
end do
diff --git a/fortran/src/read_fits4d.f90 b/fortran/src/read_fits4d.f90
index 44464c5..9d11889 100644
--- a/fortran/src/read_fits4d.f90
+++ b/fortran/src/read_fits4d.f90
@@ -59,7 +59,7 @@ subroutine read_fits4d(file,xvin,yvin,profile,lsprof,do_ls,ok)
if (fits4dfile.ne.file) then
if ((ipt%verbose.ge.1).and.(myid.eq.0)) write(*,*) 'Reading Fits 4D file ',trim(file)
- ! Get an unused Logical Unit Number to use to open the FITS file.
+ ! Get an unused Logical Unit Number to use to open the FITS file.r
call ftgiou(unit,fstatus)
! Open the FITS file previously created by WRITEIMAGE
call ftopen(unit,file,0,blocksize,fstatus)
diff --git a/idlpro/atminput_fill.pro b/idlpro/atminput_fill.pro
index 31b24ac..c3c0a38 100644
--- a/idlpro/atminput_fill.pro
+++ b/idlpro/atminput_fill.pro
@@ -6,15 +6,16 @@ function atminput_fill,xv,yv
pnam=def_parname(tag_names(ipt.atm.par))
for ia=0,ipt.ncomp-1 do begin
for i=0,atminput[ia].npar-1 do begin
- if ((xv gt atminput[ia].nx) or (yv gt atminput[ia].ny)) then begin
+ if ((xv ge atminput[ia].nx) or (yv ge atminput[ia].ny)) then begin
print,'ATMINPUT_FILL: Cannot find x=',xv,'y=',yv,' in ',ipt.atm[ia].atm_input
endif else begin
it=(where(strtrim(pnam.ipt) eq strtrim(atminput[ia].par[i])))[0]
if it eq -1 then begin
- message,'Could not find parameter from ATMINPUT: '+atminput[ia].par[i]
+ if ipt.verbose ge 1 then $
+ message,/cont,'Could not find parameter from ATMINPUT: '+atminput[ia].par[i]
endif else begin
+ retipt.atm[ia].par.(it)=atminput[ia].data[i,xv,yv]
endelse
- retipt.atm[ia].par.(it)=atminput[ia].data[i,xv,yv]
endelse
endfor
endfor
diff --git a/idlpro/helix.pro b/idlpro/helix.pro
index fd774a1..a54eeee 100644
--- a/idlpro/helix.pro
+++ b/idlpro/helix.pro
@@ -336,6 +336,7 @@ pro helix,ifile,ipt=input_file,list=plist,savall=savall, $
endif
;get initial condition / sythesis info from ATMINPUT keyword.
+; ipt=atminput_fill(xp-ipt.x[0],yp-ipt.y[0])
ipt=atminput_fill(xp,yp)
;use initial value as provided with keyword
diff --git a/idlpro/read_atminput.pro b/idlpro/read_atminput.pro
index cb09f6e..0999621 100644
--- a/idlpro/read_atminput.pro
+++ b/idlpro/read_atminput.pro
@@ -13,10 +13,17 @@ pro read_atminput
print,format='(a,i2,'': '',a)','Reading ATM_INPUT fits file for COMP ',ia, $
aistr[ia]
- data=read_fits_dlm(aistr[ia],0,header)
+ data=read_fits_dlm(aistr[ia],0,header,status=status)
+ if status ne 0 then message,'Could not find ATM_INPUT file: '+aistr[ia]
+ ;if file is an atm_ file form HeLIx+
+ ;then do the transpose
+ if strpos(sxpar(header,'CODE'),'HELIX+') eq 0 then begin
+ data=transpose(data,[2,0,1])
+ endif
naxes=size(data,/dim)
dummy=execute('data'+n2s(ia)+'=temporary(data)')
+
if ipt.verbose ge 2 then print,'npar=',naxes[0],' nx=',naxes[1],' ny=',naxes[2]
if (ia eq 0) then begin
nx=naxes[1]
--
GitLab