#
# Makefile designed for use with G77, the GNU FORTRAN 77 compiler.
#
# Copyright (C) 2005-2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
#

#
# Fortran 77-specific declarations:
#
F77     = g77
F77OPTS = -Wall -pedantic

#
# AWK-specific declarations:
#
AWK     = awk
METAPOST= mpost
TEX     = tex
DVIPS   = dvips

#
# Declarations specific for the sample execution procedure.
#
ALPHASTART=-10.0
ALPHASTOP=5.0
BETA=10.0
ETA=3.0
LENGTH=3.0
KAPPA=1.4

all: mopagras zgaussjdrv ccubsolvdrv
	@make sample
	@make linear

mopagras: mopagras.o mopagrat.o ccubsolv.o zgaussj.o
	$(F77) $(F77OPTS) -o mopagras mopagras.o mopagrat.o ccubsolv.o zgaussj.o

mopagras.o: mopagras.f
	$(F77) $(F77OPTS) -c mopagras.f

mopagrat.o: mopagrat.f
	$(F77) $(F77OPTS) -c mopagrat.f

ccubsolv.o: ccubsolv.f
	$(F77) $(F77OPTS) -c ccubsolv.f

zgaussj.o: zgaussj.f
	$(F77) $(F77OPTS) -c zgaussj.f

zgaussjdrv: zgaussjdrv.o zgaussj.o
	$(F77) $(F77OPTS) -o zgaussjdrv zgaussjdrv.o zgaussj.o

zgaussjdrv.o: zgaussjdrv.f
	$(F77) $(F77OPTS) -c zgaussjdrv.f

ccubsolvdrv: ccubsolvdrv.o ccubsolv.o
	$(F77) $(F77OPTS) -o ccubsolvdrv ccubsolvdrv.o ccubsolv.o

ccubsolvdrv.o: ccubsolvdrv.f
	$(F77) $(F77OPTS) -c ccubsolvdrv.f

sample:
	@make sample-data
	@make sample-data-postprocessing
	@make sample-generate-data-for-nonmodulated-medium
	@make sample-generate-data-for-ideally-phasematched-medium
	@make sample-graph

#
#	@./mopagras --verbose -o sample.dat \
#	   --signal_lcp 1.0 --idler_lcp 0.0 \
#	   --alpha  -15.0.0:15.0:900 \
#	   --beta   10.0:0.0:1 \
#	   --gamma  0.0:0.0:1 \
#	   --delta  0.0:0.0:1 \
#	   --eta    3.05:1.0:1 \
#	   --kappa  0.8:5.0:1 \
#	   --length 3.0:1.0:1 \
#	   --zeta   3.0:1.0:1
#
sample-data:
	@echo "-----------------------------------------------------------"
	@echo "sample :: Generating raw data for parametric conversion."
	@echo "-----------------------------------------------------------"
	@./mopagras --verbose \
	   --outputfile_lcp sample.lcp.dat \
	   --outputfile_rcp sample.rcp.dat \
	   --signal_lcp 1.0 --idler_lcp 0.0 \
	   --signal_rcp 1.0 --idler_rcp 0.0 \
	   --alpha  $(ALPHASTART):$(ALPHASTOP):900 \
	   --beta   $(BETA):$(BETA):1 \
	   --gamma  0.0:0.0:1 \
	   --delta  0.0:0.0:1 \
	   --eta    $(ETA):$(ETA):1 \
	   --kappa  $(KAPPA):$(KAPPA):1 \
	   --length $(LENGTH):$(LENGTH):1 \
	   --zeta   $(LENGTH):$(LENGTH):1

sample-data-postprocessing:
	@echo "-----------------------------------------------------------"
	@echo "sample :: Numerical postprocessing of simulated data."
	@echo "-----------------------------------------------------------"
	@$(AWK) '\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	{\
	   printf("%-2.8f %-2.8f\n",$$1,10.0*log10($$2));\
	}' sample.lcp.dat>sample.lcp.a1f.dat
	@$(AWK) '\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	{\
	   printf("%-2.8f %-2.8f\n",$$1,10.0*log10($$3));\
	}' sample.lcp.dat>sample.lcp.a1b.dat

sample-generate-data-for-nonmodulated-medium:
	@echo "-----------------------------------------------------------"
	@echo "sample :: Generating data corresponding to nonmodulated media"
	@echo "-----------------------------------------------------------"
	@echo 1|$(AWK) '\
	   function cosh(x  ,tmp) {\
	      tmp=((exp(x)+exp(-x))/2.0);\
	      return tmp;\
	   }\
	   function sinh(x  ,tmp) {\
	      tmp=((exp(x)-exp(-x))/2.0);\
	      return tmp;\
	   }\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	BEGIN {\
	   alphastart='$(ALPHASTART)';\
	   alphastop='$(ALPHASTOP)';\
	   beta='$(BETA)';\
	   eta='$(ETA)';\
	   ll=3.0;\
	   for(k=1;k<=2;k++) {\
	      alpha=alphastart+(k-1)*(alphastop-alphastart);\
	      tmp=eta*eta-beta*beta/4.0;\
	      if (tmp>0.0) {\
	         kk=sqrt(tmp);\
	         tmp=cosh(kk*ll)*cosh(kk*ll);\
	         tmp=tmp+(beta*beta/(4.0*kk*kk))*sinh(kk*ll)*sinh(kk*ll);\
	      } else {\
	         kk=sqrt(-tmp);\
	         tmp=cos(kk*ll)*cos(kk*ll);\
	         tmp=tmp+(beta*beta/(4.0*kk*kk))*sin(kk*ll)*sin(kk*ll);\
	      }\
	      tmp=10.0*log10(tmp);\
	      printf("%-2.8f %-2.8f\n",alpha,tmp);\
	   }\
	}\
	END {}' sample.lcp.dat>sample.lcp.a1f.ref.dat

sample-generate-data-for-ideally-phasematched-medium:
	@echo "-----------------------------------------------------------"
	@echo "sample :: Generating data corresponding to ideal phase matching"
	@echo "-----------------------------------------------------------"
	@echo 1|$(AWK) '\
	   function cosh(x  ,tmp) {\
	      tmp=((exp(x)+exp(-x))/2.0);\
	      return tmp;\
	   }\
	   function sinh(x  ,tmp) {\
	      tmp=((exp(x)-exp(-x))/2.0);\
	      return tmp;\
	   }\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	BEGIN {\
	   alphastart=-15.0;\
	   alphastop=15.0;\
	   beta='$(BETA)';\
	   eta='$(ETA)';\
	   ll=3.0;\
	   for(k=1;k<=2;k++) {\
	      alpha=alphastart+(k-1)*(alphastop-alphastart);\
	      tmp=eta*eta-beta*beta/4.0;\
	      if (tmp>0.0) {\
	         kk=sqrt(tmp);\
	         tmp=cosh(kk*ll)*cosh(kk*ll);\
	         tmp=tmp+(beta*beta/(4.0*kk*kk))*sinh(kk*ll)*sinh(kk*ll);\
	      } else {\
	         kk=sqrt(-tmp);\
	         tmp=cos(kk*ll)*cos(kk*ll);\
	         tmp=tmp+(beta*beta/(4.0*kk*kk))*sin(kk*ll)*sin(kk*ll);\
	      }\
	      tmp=10.0*log10(tmp);\
	      printf("%-2.8f %-2.8f\n",alpha,tmp);\
	   }\
	}\
	END {}' sample.lcp.dat>sample.lcp.a1f.pm.dat

sample-graph:
	@echo "-----------------------------------------------------------"
	@echo "sample :: Creating graph of signal conversion efficciency"
	@echo "-----------------------------------------------------------"
	@echo \
	   "input graph; beginfig(1); w:=80mm; h:=56mm;"\
	   "draw begingraph(w,h);"\
	   "pickup pencircle scaled .3pt;"\
	   "setrange(whatever,whatever,whatever,whatever);"\
	   "pickup pencircle scaled .5pt;"\
	   "gdraw \"sample.lcp.a1f.dat\";"\
	   "gdraw \"sample.lcp.a1f.ref.dat\" dashed evenly;"\
	   "gdraw \"sample.lcp.a1b.dat\" dashed evenly;"\
	   "pickup pencircle scaled .3pt;"\
	   "glabel.bot(btex\hbox to 7mm{\hfil} $$ \alpha$$ etex,OUT);"\
	   "glabel.lft(btex$$\
	      10\log_{10}(|a^{{\rm f}\pm}_{1}(L)/a^{{\rm f}\pm}_{1}(0)|^2)\
	      $$ etex rotated 90,OUT);"\
	   "endgraph; endfig; end" > sample.lcp.mp
	@$(METAPOST) sample.lcp.mp
	@$(TEX) -jobname=sample.lcp '\input epsf\nopagenumbers\
	   \centerline{\epsfxsize=120mm\epsfbox{sample.lcp.1}}\bye'
	@$(DVIPS) -D1200 -E sample.lcp.dvi -o sample.lcp.eps

linear:
	@make linear-data
	@make linear-analytic-solution
	@make linear-data-postprocessing
	@make linear-graph
	@make linear-log-graph

linear-data:
	@echo "-----------------------------------------------------------"
	@echo "linear :: Generating raw data for parametric conversion."
	@echo "-----------------------------------------------------------"
	@./mopagras --verbose \
	   --outputfile_lcp linear.lcp.dat \
	   --outputfile_rcp linear.rcp.dat \
	   --signal_lcp 1.0 --idler_lcp 0.0 \
	   --signal_rcp 1.0 --idler_rcp 0.0 \
	   --alpha  -10.0:10.0:900 \
	   --beta   15.0:10.0:1 \
	   --gamma  0.0:0.0:1 \
	   --delta  0.0:0.0:1 \
	   --eta    0.001:0.001:1 \
	   --kappa  3.0:3.0:1 \
	   --length 2.0:2.0:1 \
	   --zeta   2.0:2.0:1

linear-analytic-solution:
	@echo "-----------------------------------------------------------"
	@echo "linear :: Generating data corresponding to analytic solution"
	@echo "-----------------------------------------------------------"
	@echo 1|$(AWK) '\
	   function cosh(x  ,tmp) {\
	      tmp=((exp(x)+exp(-x))/2.0);\
	      return tmp;\
	   }\
	   function sinh(x  ,tmp) {\
	      tmp=((exp(x)-exp(-x))/2.0);\
	      return tmp;\
	   }\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	BEGIN {\
	   alphastart=-10.0;\
	   alphastop=10.0;\
	   kappa=3.0;\
	   ll=2.0;\
	   nn=600;\
	   for(k=1;k<=nn;k++) {\
	      alpha=alphastart+((k-1)*(alphastop-alphastart))/(nn-1);\
	      tmp=kappa*kappa-alpha*alpha;\
	      if (tmp>0.0) {\
	         xi=sqrt(tmp);\
	         tmp=xi*xi*cosh(xi*ll)*cosh(xi*ll);\
	         tmp=tmp+alpha*alpha*sinh(xi*ll)*sinh(xi*ll);\
	         tmp=xi*xi/tmp;\
	      } else {\
	         xi=sqrt(-tmp);\
	         tmp=xi*xi*cos(xi*ll)*cos(xi*ll);\
	         tmp=tmp+alpha*alpha*sin(xi*ll)*sin(xi*ll);\
	         tmp=xi*xi/tmp;\
	      }\
	      printf("%-2.8f %-2.8f\n",alpha,tmp);\
	   }\
	}\
	END {}' linear.lcp.dat>linear.lcp.tt.ref.dat

linear-data-postprocessing:
	@echo "-----------------------------------------------------------"
	@echo "linear :: Numerical postprocessing of simulated data."
	@echo "-----------------------------------------------------------"
	@$(AWK) '{printf("%-2.8f %-2.8f\n",$$1,$$2)}' linear.lcp.dat>linear.lcp.tt.dat
	@$(AWK) '{printf("%-2.8f %-2.8f\n",$$1,$$3)}' linear.lcp.dat>linear.lcp.rr.dat
	@$(AWK) '\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	{\
	   printf("%-2.8f %-2.8f\n",$$1,10.0*log10($$2));\
	}' linear.lcp.dat>linear.lcp.ttdb.dat
	@$(AWK) '\
	   function log10(x) {\
	      return log(x)/log(10.0);\
	   }\
	{\
	   printf("%-2.8f %-2.8f\n",$$1,10.0*log10($$2));\
	}' linear.lcp.tt.ref.dat>linear.lcp.ttdb.ref.dat

linear-graph:
	@echo "-----------------------------------------------------------"
	@echo "linear :: Creating graph of linear transmission of grating"
	@echo "-----------------------------------------------------------"
	@echo \
	   "input graph; beginfig(1); w:=80mm; h:=56mm;"\
	   "draw begingraph(w,h);"\
	   "pickup pencircle scaled .3pt;"\
	   "setrange(whatever,whatever,whatever,whatever);"\
	   "pickup pencircle scaled .5pt;"\
	   "gdraw \"linear.lcp.tt.dat\";"\
	   "gdraw \"linear.lcp.tt.ref.dat\" dashed evenly;"\
	   "gdraw \"linear.lcp.rr.dat\" dashed evenly withcolor white;"\
	   "pickup pencircle scaled .3pt;"\
	   "glabel.bot(btex\hbox to 7mm{\hfil} $$ \alpha$$ etex,OUT);"\
	   "glabel.lft(btex$$\
	      10\log_{10}(|a^{{\rm f}\pm}_{1}(L)/a^{{\rm f}\pm}_{1}(0)|^2)\
	      $$ etex rotated 90,OUT);"\
	   "endgraph; endfig; end" > linear.lcp.mp
	@$(METAPOST) linear.lcp.mp
	@$(TEX) -jobname=linear.lcp '\input epsf\nopagenumbers\
	   \centerline{\epsfxsize=120mm\epsfbox{linear.lcp.1}}\bye'
	@$(DVIPS) -D1200 -E linear.lcp.dvi -o linear.lcp.eps

linear-log-graph:
	@echo "-----------------------------------------------------------"
	@echo "linear :: Creating graph of linear transmission of grating"
	@echo "-----------------------------------------------------------"
	@echo \
	   "input graph; beginfig(1); w:=80mm; h:=56mm;"\
	   "draw begingraph(w,h);"\
	   "pickup pencircle scaled .3pt;"\
	   "setrange(whatever,whatever,whatever,whatever);"\
	   "pickup pencircle scaled .5pt;"\
	   "gdraw \"linear.lcp.ttdb.dat\";"\
	   "gdraw \"linear.lcp.ttdb.ref.dat\" dashed evenly;"\
	   "pickup pencircle scaled .3pt;"\
	   "glabel.bot(btex\hbox to 7mm{\hfil} $$ \alpha$$ etex,OUT);"\
	   "glabel.lft(btex$$\
	      T=|a^{{\rm f}\pm}_{1}(L)/a^{{\rm f}\pm}_{1}(0)|^2\
	      $$ etex rotated 90,OUT);"\
	   "endgraph; endfig; end" > linear.lcp.db.mp
	@$(METAPOST) linear.lcp.db.mp
	@$(TEX) -jobname=linear.lcp.db '\input epsf\nopagenumbers\
	   \centerline{\epsfxsize=120mm\epsfbox{linear.lcp.db.1}}\bye'
	@$(DVIPS) -D1200 -E linear.lcp.db.dvi -o linear.lcp.db.eps

clean:
	-rm -f *~ *.o *.exe *.dat *.stackdump mopagras *.mp *.mpx *.tex
	-rm -f *.1 *.dvi *.log *.eps zgaussjdrv cubctest ccubsolvdrv

archive:
	make -ik clean
	tar --gzip --directory=../ -cf mopagras.tgz mopagras